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INTRODUCTION 


Practical,  mass-produced  rocket  propellants  usually  contain  trace 
amounts  of  impurities  such  as  potassium,  sodium,  or  calcium.  Even  with 
very  low  concentrations  (~  100  ppm),  these  low  ionization  potential 
species  are  sufficiently  ionized  to  contribute  a  substantial  number  of 
free  electrons  to  the  exhaust.  In  the  presence  of  these  impurities  or 
additives,  the  rocket  motor  exhaust  flow  becomes  a  weakly  ionized 
plasma.  The  exhaust  plume  becomes  electrically  active  and  affects  the 
propagation  of  those  electromagnetic  (EM)  waves  with  wavelengths  in  the 
microwave  r«nge  or  longer.  The  exhaust  plume  can  act  to  attenuate  and 
scatter  incident  waves  which  are  intended  to  carry  information  to  or 
from  the  missile.  The  exhaust  plume  reflects  and  scatters  waves  which 
contribute  to  a  radar  cross  section  in  addition  to  the  contribution  of 
the  missile  body.  This  class  of  interactions  is  the  exhaust  plume 
radar  frequency  (RF)  interference  problem. 

Plume  RF  interference  has  an  extensive  history  because  it  is  a 
potential  problem  for  nearly  every  missile  system.  Past  problems  have 
been  resolved  by  a  judicious  combination  of  experiments  and  analysis. 
A  comprehensive  summary  of  this  work  has  been  published.1-3  This 
earlier  work  identifies  the  importance  of  two  major  components  of  the 
problem:  (1)  the  exhaust  plume  flowfield  and  composition,  and  (2)  the 
interaction  of  the  EM  wave  with  the  flow.  There  has  recently  been  a 
step  forward  in  our  capability  to  deal  with  the  first  issue.  The  first 
version  of  the  JANNAF  Standard  Plume  Flowfield  (SPF/1)  Model  has  been 
released  to  industry.1*  This  code  allows  an  improved  description  of  the 
exhaust  plume  that  can  account  for  both  the  inertially  dominated  region 
of  the  plume  (initial  expansion)  and  the  turbulent  mixing  with  the 


Wictor,  A.C. ,  "Plume-Signal  Interference,  Part  1.  Radar  Attenuation," 
Naval  Weapons  Center,  China  Lake,  Calif.,  June  1975.  (NWC  TP  5319, 

publication  UNCLASSIFIED.) 

i - .  "Plume-Signal  Interference,  Part  2.  Plume-Induced  Noise," 

Naval  Weapons  Center,  China  Lake,  Calif.,  May  1972.  (NWC  TP  5319, 

publication  formerly  CONFIDENTIAL,  declassified  5  April  1977.) 

3JANNAF  Handbc.  >k.  Rocket  Exhaust  Plume  Technology,  Chapter  4.  Plume 
Electromagnet!.  Interactions,  April  1977.  (CPIA  Pub.  263,  publication 

unclassified;) 

'’Dash,  S.M.  and  Pergament,  H.S.,  "The  JANNAF  Standard  Plume  Flowfield 
Model  (SPF),"  Aeronautical  Research  Associates  of  Princeton,  Inc., 
N.J.,  April  1981.  (A.R.A.P.  Report  No.  448,  publication  UNCLASSIFIED. ) 
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ambient  air.  Very  general  chemical  kinetics  are  included.*  Moreover, 
there  is  a  potential  to  provide  some  initial  predictions  of  turbulence 
quantities  of  importance  to  scattering  of  RF  waves  by  extension  of  the 
SPF  to  include  equations  for  the  evolution  of  the  second-order 
correlations  of  passive  scalar  quantities  (i.e.f  the  g-equation) 

This  report  describes  a  step  forward  in  dealing  with  the  second 
category  of  the  interaction  of  the  EM  wave  with  the  flow.  It  describes 
the  Naval  Weapons  Center  Plume  Radar  Frequency  Interference  Code 
(PRFIC)  which  is  intended  to  unify  the  calculations  performed  by 
previous,  separate  codes  and  to  utilize  the  flowfleld  as  predicted  by 
the  JANNAF  SPF.  The  combination  of  these  two  codes  offers  the 
potential  for  improved  predictions  of  RF  plume  interference. 

This  report  contains  a  physical  description  of  plume  RF  phenomena 
and  includes  examples  of  both  the  flowfleld  calculation  and  EM 
interaction  predictions.  The  PRFIC  user's  manual  is  provided  in 
Appendix  A,  and  the  software  documentation,  in  Appendix  B. 


PROPAGATION  OF  AN  EM  WAVE  IN  A  WEAKLY  IONIZED  PLASMA 


Propagation  of  an  electromagnetic  wave  through  an  exhaust  plume  is 
modeled  by  a  plane,  periodic  wave  of  the  form6 

E  ~  elkx  e-i“t  (1) 


The  wavenumber  k  is  complex  in  general.  We  express  it  as  k  »  S  +  ia, 
where  3  is  a  phase  constant,  and  u,  the  imaginary  part,  is  an 
attenuation  constant.  If  this  form  of  the  wave  is  substituted  into 
Maxwell's  equations,  the  resulting  dispersion  relation  is  obtained. 


*In  the  laminar  limit,  i.e.,  there  is  no  effect  of  turbulence  on 
reaction  rates  other  than  the  influence  of  the  mean  flow. 

sKhalil,  E.E.,  Scalding,  O.B.,  and  Whitelaw,  J.H. ,  "The  Calculation  of 
Local  Flow  Properties  in  Two-Dimensional  Furnaces,”  Int.  J.  Heat  Mass 
Transfer .  Vol.  18,  1975,  pp.  775-791. 

6Mitchner,  M.  and  Kruger,  C.  H. ,  Partially  Ionized  Gases.  John  Wiley 
and  Sons,  N.Y. ,  1973,  pp.  156-61. 
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(2) 


where 


<*>p  *  (nee2/eome)1/2  (?) 

is  the  plasma  frequency  (Vp  =  u)p/2ir  *  8.97  n^/2  Hz,  ne  in  n-3).  ven  is 
the  electron-neutral  collision  frequency 


cn 


N 

ven  L  xi^i 
i-1 


(4) 


where  ve  is  the  electron  thermal  velocity  ( 3kT/me ) 1 1 2 ,  n  is  the  total 
particle  number  density  p/kT,  and  x^  is  the  mole  fraction  of  species  i. 

is  the  collision  cross  section,  which  can  depend  on  the  colliding 
molecule  and  temperature  (or  equivalently,  electron  velocity).  The 
summation  is  over  the  total  number  of  species  i  =  1....N.  For  the 

calculation  given  here,  we  used  «  const  »  1.5xl0~14  m2.  It  is 
useful  to  note  that  ven  -  p  T-1/2. 


ATTENUATION  AND  PHASE 


The  plasma  frequency  is  the  natural  frequency  of  oscillation  of 
electrons  if  displaced  from  their  equilibrium  positions  in  the  force 
field  of  the  heavier  ions.  In  the  absence  of  collisions,  the 
attenuation  of  a  wave  depends  discontinuously  on  whether  this  natural 
frequency  is  greater  or  less  than  the  wave  frequency.  For  the  special 
case  of  no  collisions,  if  u)p  >  w,  the  electrons  can  respond 
sufficiently  rapidly  to  neutralize  the  incident  wave.  k  is  purely 
imaginary, 

(5) 
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and  no  wave  of  the  form  of  Ea .  (1)  can  propagate  Into  the  plasma.  This 
is  the  "overdense"  case  in  which  the  electrons  can  adjust  rapidly  to  a 
disturbance  to  reflect  the  incident  energy.  Again,  for  the  special 
case  of  few  collisions,  the  electrons  cannot  respond  quickly  to  cancel 
the  incident  wave  if  ujp  <  w;  then  k  is  purely  real. 


and  there  is  no  attenuation  of  the  wave.  This  is  the  "underdense" 
case.  With  finite  collision  frequency,  the  electron  oscillations  are 
damped  and  the  discontinuous  nature  of  the  propagation  at  j  *  u)p  is 
smoothed.  There  is  propagation  and  attenuation  throughout  the  entire 
range  of  j. 

The  attenuation  and  phase  constants  are  given  explicitly  if  we  let 

-  *.  +  1  Ki  <7> 
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Normalized  values  of  the  attenuation  and  phase  constants  (uc/m  and 
6c/j)  are  given  in  Figures  l  and  2,  respectively,  for  a  range  of  u>p/u> 
with  collision  frequency  (ven/w)  as  a  parameter.  As  the  collision 

frequency  increases,  the  transition  through  the  critical  electron 
number  density  at  which  u>  *  ujp  becomes  less  abrupt,  and  the  notion  of 
overdense  and  underdense  becomes  less  meaningful.  It  is  important  to 
note  that  the  maximum  attenuation  at  a  given  frequency  ratio  (oop/M) 
occurs  for  ven/ti)  of  order  unity. 

Magnitudes  of  the  plasma  frequency  (vp  *  jp/2-n),  the  collision 

frequency  ven,  and  the  attenuation  constant  u(in  m-^)  are  shown  in 
contour  plots  for  a  typical  exhaust  plume  in  Figures  3,  A,  and  5. 

These  predictions  are  for  76  percent  solids,  12  percent  aluminum 

propellent,  in  a  20,000  lbf  thrust  motor  flying  at  a  sonic  speed  at  2 
km  altitude.  This  plume  is  underdense  to  all  radar  frequencies  greater 
than  _.n  GHz.  However,  the  ratio  of  collision  frequency  to  wave 
frequency  is  large  (ven/ui  >  1).  Figure  5  suggests  that  diagonal 

attenuation  through  this  plume  will  be  severe. 

It  is  useful  to  consider  the  limiting  behavior  of  the  attenuation 

and  phase  for  the  extremes  jf  electron  number  density  and  collision 

frequency.  These  limits  are  defined  in  Table  1.  Low  altitude,  solid 
propellant  motor  exhaust  plumes  correspond  to  ven/uj  >>  1  and  -Vp/w  <  1. 
Liquid  propellant  exhaust  plumes  usually  have  lower  impurity  levels 
such  that  !j)p/<i)<<l.  The  attenuation  and  phase  shift  of  wave  propagation 
through  the  plumes  of  missiles  therefore  depend  strongly  on  the  class 
of  propellant  and  altitude. 

TABLE  1.  Limiting  Forms  for  the  Attenuation 
and  Phase  Constants 


attenuation,  a: 


phase  ,  £  : 


Normalized  Attenuation  Constant,  ac/<u 
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FIGURE  1.  Attenuation  Constant  in  an  Ionized 
Gas  with  Col  1 i sions  . 
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The  definition  of  the  attenuation  constant  in  Eq .  (1)  applies  to 
the  reduction  of  the  electric  field.  In  RF  interference  problems,  the 
usual  measurement  is  of  the  power,  which  depends  on  E2.  The 
attenuation  used  in  the  remainder  of  this  report  is  therefore  2a,  and 
the  attenuation  is  given  in  dB.* 

Attenuation  and  the  phase  change  along  an  optical  path  through  the 
plume  involve  an  integral  of  a  and  S  along  the  optical  path  3, 


* 


p"2ux 


P  P 

—  .  dB  ■  10  loJtjQ  — •  *  2()ax  log^  e  «  8.659ax 
Po  Po 
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attenuat Ion 


(12: 


phase 


(13) 


In  general,  this  optical  path  is  curved.  The  propagation  is  through  a 
medium  with  variable  refractive  index,  and  the  plume  acts  as  a  lens  at 
these  wavelengths.  In  the  geometric  optics  limit,  the  ray  path  is 
described  by/ 


_d_ 

ds 


(nt  ) 


V  n 


(14) 


t  is  the  unit  vector  tangent  to  the  ray  path,  and  s  is  the  arc  length 
along  the  path.  This  path  is  curved  if  there  is  a  component  of  Vn 
normal  to  the  direction  of  propagation.  Eq .  (14)  is  equivalent  to 
three  ordinary  differential  equations  for  the  coordinate  of  the  ray 
path  x(s),  y(s),  z(s).  In  the  case  of  an  axisymmetric  exhaust  plume, 
there  are  only  two  components  of  Vn:  in  the  radial  and  axial 
directions.  The  complete  system  to  be  solved  is  (see  Apnendix  C  for 
the  derivation  from  Eq .  (14)). 


d^x 

ds^ 


cos9 


n 


sind 


n 


3n  dx  / 

-  T~  COSU 

3r  ds  y 

3n  dy  / 

3r  '  ds(C°Sa 


dx 

ds 


dx 

ds 


+  sinO  d-y-  3-°-  +  ^ 

3r  ds  dr  ds  dz/ 


3n 

— -  +  sind 
3r 


d_y  3n  dz  3n 
ds  3r  ds  3z 


(15) 


(16) 


d2z  1  3n  dz  /  dx  3n  .  dy  3n  dz  ?n\ 

- — j  »  -  —  -  —  cosy  —  —  +  sine  —  —  +  --  — 

ds‘  n  3z  ds  1  ds  3r  ds  3r  ds  3z / 


(17) 


In  general,  we  know  the  initial  point  x(0),  y(0),  z(0),  and  the  initial 
direction  of  the  ray  dx/ds,  dy/ds,  dz/ds  at  s  »  0.  The  solution 
provides  the  final  location  and  direction  after  traversing  the  region 


^Marchand,  E.  W.  ,  Gradient  Index  Optics,  Academic  Press,  N.Y.,  197H, 

P.  108.  "  ' 
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of  variable  Index  of  refraction.  The  index  of  refraction  at  any  point 
is 


n  *  cg/ui 


(18) 


In  free  space,  n  =■  1 ,  and  g  *  2ir/X.  In  the  limit  of  small  gradients  in 
n,  or  for  paths  along  which  t  •  Vn  ■  0,  the  ray  trajectory  is  a 
straight  line-of-sight  (LOS). 

In  our  numerical  evaluation  of  this  path.  the  differential 
equations  are  solved  by  finite  differences. 


As2  (19) 

s 


x(s  +  As) 


x(s)  +  ~ 
ds 


,  A  1  d2x 

*S  +  2  dP- 


y(s  +  As) 


y(s) 


s 


As  + 


1  l£z 

2  ds2 


As2 


(20) 


z(s  +  As) 


z(s) 


As  + 


A  Al?. 

2  ds2 


As2 

s 


(21) 


A  solution  proceeds  from  the  known  initial  location  x(0),  y(0),  z(0) 
with  known  direction  [dx(0)/ds],  [dy(0)/ds],  [dz(0)/dx]  in  a  stepwise 
manner  to  the  final  point  x(s),  y(s),  z(s).  However,  this  final  point 
is  known  in  advance  (the  receiver  location),  but  the  initial  direction 
of  the  refracted  ray  required  to  intercept  the  final  point  is  unknown. 
This  is  a  two-point  boundary  value  problem  in  which  the  initial 
direction  which  allows  the  trajectory  to  intercept  the  end  point  is  to 
be  found.  The  numerical  solution  proceeds  by  choosing  an  initial  value 
for  the  direction  at  the  transmitter  (the  direction  of  the 
line-of-sight),  and  then  iterating  until  the  final  location  is  within  a 
prescribed  distance  of  the  receiver  location.  Each  integration 
requires  a  new  guess  of  the  direction,  which  we  have  taken  to  be 

(cos  ai)(n+1)  -  (cos  at)(n)  +  RC  •  e^s  (22) 
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where  is  the  direction  cosine  with  respect  to  the  i^  direction,  €< 
is  the  component  of  the  error  vector  from  the  end  point  of  the  nc“ 
integration,  and  s  is  the  ray  path  length  to  that  point.  RC  is  an 
overrelaxation  factor  (0  <  RC  <  1)_.  Thu  end  of  each  integration  is 
taken  to  be  the  point  at  which  |e|,  the  error  between  the  n1-*1 
integration  and  the  actual  end  point,  is  a  minimum  (i.e.,  the  error  is 
computed  at  each  As,  and  the  integration  stops  when  |c|  begins  to 
increase).  When  the  actual  path  is  obtained,  |e|  <  e^,  and  the 
attenuation  and  phase  are  evaluated  along  the  curved  path,  s,  according 
to  Eqs.  ( 12)  and  (13). 


SCATTERING  BY  TURBULENCE 


The  SPF  accounts  for  turbulence  in  the  exhaust  plume  flow.  The 
primary  emphasis  is  on  mixing  and  specification  of  the  mean  flowfield 
properties,  but  predictions  of  turbulence  quantities  are  also 
performed.  For  example,  a  dynamical  equation  for  the  evolution  of  the 
turbulence  kinetic  energy  (k  *  1/2  q ' ^ )  is  solved  in  the  two-equation 
turbulence  models.  Improved  descriptions  involving  the  g-equation 
formulation  of  the  mean  square  fluctuations  of  passive  scalar 
quantities  are  currently  included  in  preliminary  form.  The  potential 
of  these  models,  and  the  even  more  detailed  turbulence  descriptions 
available  in  higher-order  closure  models,  are  just  beginning  to  be 
realized  in  exhaust  plume  applications. 


DIFFERENTIAL  CROSS  SECTION 

In  view  of  the  potential  availability  of  plume  solutions  Including 
a  turbulence  description,  PRFIC  includes  a  formulation  for  scattering 
from  turbulent  fluctuations  of  the  index  of  refraction.  We  have  used 
the  single  scattering  (Born)  approximation,  but  with  an  added 
correction  to  account  foi  attenuation  by  the  mean  flow  of  the  incident 
and  scattered  wave. 

*It  is  perhaps  useful  to  point  out  the  distinction  between  scattering 
and  diffraction.  Scattering  is  the  more  general  terra,  and  usually 
refers  to  the  distribution  of  radiation  at  distances  from  the  volume 
of  nonuniform  refractive  index  that  are  large  compared  to  the 
wavelength  and  to  the  size  of  the  nonuniform  volumes.  The  effect  of 
turbulence  on  RF  waves,  in  which  the  size  of  the  most  effective 
nonuniformities  is  comparable  to  the  wavelength  and  is  small  compared 
to  the  usual  distances  to  the  transmitter  and  receiver,  is  usually 
termed  scattering.  In  contrast,  the  effect  of  the  mean  flow,  for 
which  the  dimensions  are  usually  large  compared  to  the  wavelength  and 
can  be  comparable  to  the  distances  to  the  transmitter  or  receiver,  is 
usually  called  diffraction.  Another  distinction  is  that  scattering  is 
a  volume  phenomenon,  while  diffraction  is  a  surface  phenomenon. 
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In  the  single  scattering  limit,  the  differential  cross  section  per 
unit  volume  of  scattering  material  isB»^ 


d2o 

dVdfl 


2  it  k>  sin2x$(S  -  k 


i) 


(23) 


$  is  the  power  spectral  density  of  the  index  of  refraction 
fluctuations,  evaluated  at  the  wavenumber  which  is  the  difference 
between  the  wavenumber  of  the  incident  and  scattered  wave. 

|k  -  k*m,  =  2k  sin  '72,  where  '  is  the  scattering  angle. 

'we  have  used  the  well-known  Kolmogorov  form  of  the  energy  spectrum 
for  which 


$  (fc  -  !c 


m) 


0.033  C* 

[l  +  4k2 A2  sin2  |  J  U/& 


(24) 


and  the  structure  constant  C2  is 


(25) 


We  use  a  length  scale  A  given  by 


A 


1  si  _L_ 

8  e  2/2 


(26) 


8Booker,  H.G. ,  and  Gordon,  W.E. ,  "A  Theory  of  Radio  Scattering  in  the 
Troposphere,"  Proceedings  of  the  I.R.E,  April  1950,  pp.  401-412. 

^Tatarski,  V.I.,  Wave  Propagation  in  a  Turbulent  Medium,  McGraw-Hill, 
New  York,  1961,  p.  68. 
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where  k  is  the  turbulence  kinetic  energy  and  e  the  dissipation  (both 
quantities  are  calculated  in  the  SPK). 

If  we  assume  that  n'  is  due  only  to  the  fluctuations  in  ,  (see 
Appendix  D),  then 


£ 7 - ^—7 

4(u>2  +  v  2) 
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(27) 


4n  2 

IF- 


1  + 


VenW 
w  / 


(28) 


with  rt  being  the  Thompson  radius  of  the  electron 


c  Afte^mc2 


2.81  x  10*15  m 


(29) 


With  these  assumptions,  the  complete  cross  section  per  unit  volume 


is 
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(30) 


The  ratio  of  received  scattered  power  to  transmitted  power  is  obtained 
from  an  integration  over  the  entire  plume  volume 
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G-p  and  GR  are  transmitter  and  receiver  gains  in  the  directions  $p  and 
$R  between  the  antenna  direction  and  the  wave  direction. 
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This  is  the  standard  formulation  for  the  power  received  from  a 
volume  of  distributed  scatterers.  We  have  modified  this  formulation  to 
include  the  attenuation  of  the  wave  along  the  path  prior  to  and 
subsequent  to  scattering  from  each  point  in  the  flow.  The  integrand  in 
Eq.  (31)  is  therefore  reduced  by  the  factor 


where  the  integrals  are  the  attenuation  along  the  optical  paths  of  the 
incident  and  scattered  waves,  respectively.  These  paths  can  be 
evaluated  along  straight  lines-of-sight  or  also  along  refracted  paths. 


ATTENUATION  DUE  TO  SCATTERING 

For  completeness,  we  note  that  scattering  contributes  to 
attenuation  by  removing  radiation  from  the  direction  of  propagation. 
It  also  reduces  attenuation  by  adding  radiation  back  into  the  direction 
of  propagation.  The  cross  section  per  unit  volume  for  removal 
(extinction)  is 
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using 


sin^x  ”  1  ~  sin20  sin^$ 


(35) 


the  integral  of  d^o/dVdW  over  all  directions  is 
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Q(p)  -  ~ 


|  (2p2  +  2p  +  1)11  -  (I  +  2p)~b/6|  -  2(p  +  1)|(1  +  2p)1/b-l] 
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Limiting  cases  are 
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The  attenuation  due  to  scattering  by  the  turbulent  fluctuations  is  then 
given  by 


n^  r{|  A3  0  (2k2A2) 

8*  3  (.0528)  -£ — 
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(39) 


This  contribution  to  'he  attenuation  along  the  line-of-sight  is 
included  in  PRF1C.  Hov^ver,  it  has  proved  to  be  small  compared  to  the 
attenuation  due  to  absorption  in  all  cases  examined  so  far. 


RADAR  CROSS  SECTION 

In  the  special  case  of  co-located  antennas,  the  scattered  power 
given  by  Eq .  (31)  reduces  to  the  turbulent  scattering  component  of  the 
radar  cross  section  of  the  plume 
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(when  the  range  R  is  large  compared  to  any  plume  dimension). 


DOPPLER  SHIFT 


SHIFT  OF  DIRECTLY  TRANSMITTED  WAVE 

Two  Doppler  shifts  of  the  transmitted  frequency  are  considered  in 
PRFIC.  One  is  the  frequency  shift  of  the  directly  transmitted  wave  due 
to  relative  motion  between  the  transmitter  and  receiver.  This  shift  is 


Av  =  v 


s 


* 


(41) 


where  s  is  a  unit  vector  from  the  transmitter  to  the  receiver.  This 
shift  occurs  for  any  condition  of  relative  motion  between  the 
transmitter  and  receiver  and  is  independent  of  the  presence  of  an 
exhaust  plume. 


TURBULENCE  BROADENED  SPECTRUM 

The  second  Doppler  shift  occurs  in  the  scattering  of  the  radiation 
by  the  turbulence  in  the  plume.  This  shift  depends  on  the  relative 
velocities  of  the  transmitter  and  receiver,  the  mean  velocity  of  the 
exhaust  pluie,  and  also  on  the  fluctuations  of  the  plume  velocity. 

Since  these  fluctuations  are  random,  they  are  described  in  a 

statistical  sense;  and,  consequently,  the  frequency  shift  due  to 

scattering  is  also  given  statistically.  In  this  first  version  of 

PRFIC,  the  velocity  fluctuations  were  assumed  to  have  a  Gaussian 
distribution  in  magnitude;  and,  consequently,  the  probability  density 
of  the  frequency  shift  due  to  scattering  is  also  Gaussian.  The 
probability  density  is 


P(Av)= 


V2^  -  •  0 


exp- 


3  2v„;^ir. 

o  2 


cAv 

Lvo 


(VT  •  s 


Vc 


t)  -  (wk  +  Vffl)(t  -  s) 


4q2"  sin2  — 
2 


(.42) 


*Terms  of  order  [(V-p  -  VR)  •  s/c|2  and  higher  are  neglected. 
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(See  Appendix  K  for  the  derivation.)  is  the  mean  square  turbulent 
velocity  fluctuation  (q^-2k,  where  k  is  the  turbulence  kinetic  energy). 
Vt,  VK.  Vm,  are  the  velocities  of  _the  transmitter,  receiver,  and 
missile,  respectively  (Figure  6).  w  is  the  local  mean  axial  velocity 
of  the  plume,  with  respect  to  a  coordinate  system  fixed  in  the  missile, 
and  is  given  directly  by  the  output  of  the  SPF. 

The  scattered  power  is  frequency  dependent,  and  the  amplitude  at 
each  frequency  v»v0+Av  is  given  by  Eq .  (31)  weighted  with  the  above 
frequency  probability  density.  Note  that  the  frequency  probability 
density  must  be  included  within  the  volume  integral  because  it  depends 
on  the  local  position  within  the  plume.  The  result  of  this  evaluation 
is  the  frequency  spectrum  of  the  scattered  power  (i.e.,  the  AM  noise 
spectrum).  The  total  AM  noise  and  total  scattering  cross  sections  are 
obtained  from 


-r—  (Av)  d(Av). 
PT 


(43) 


APPLICATIONS  OF  PRFIC 


This  section  contains  summaries  of  some  sample  calculations  with 
PRFIC  to  illustrate  its  capabilities.  As  with  any  code  that  makes 
predictions  of  physical  quantities,  the  ultimate  test  is  to  compare  the 
predictions  with  measurements.  However,  with  RF  interference,  there  is 
the  equall)  important  issue  of  the  correct  flowfield.  Comparison 
between  measured  and  predicted  RF  interference  quantities  therefore 
requires  the  additional  effort  of  making  accurate  flowfield 
predictions.  That  was  not  part  of  the  contractual  effort  summarized 
here  and  remains  to  be  done.  The  emphasis  in  this  initial  effort  was 
to  verify  that  PRFIC  computes  correctly  those  quantities  which  it  is 
supposed  to  provide,  given  a  plume  flowfield  description. 


TRANSVERSE  ATTENUATION 

We  include  one  comparison  of  predictions  and  experiment  with  the 
objective  of  showing  that  PRFIC  does  provide  correct  attenuation 
predictions.  An  extensive  set  of  transverse  attenuation  measurements 
was  reported  in  Reference  1.  We  have  made  predictions  for  cases  21-23 
of  this  set:  a  small  static,  solid  propellant  motor  (84  percent 
solids,  of  which  18  percent  is  aluminum). 


22 


NWC  TP  6336 


Receiver 


FIGURE  6.  Coordinate  System  for  PRFIC. 
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The  full  capability  of  the  SPF  was  utilized.  We  Input  a  user 
defined  chemical  system  with  species  and  reactions  that  included  all 
the  major  afterburning  and  charged  species  and  ail  the  afterburning  and 
charged  species  reactions.  The  impurity  levels  were  80  ppm  potassium 
and  126  ppra  sodium.  Both  equilibrium  and  finite  rate  chemistry  during 
the  expansion  from  the  chamber  were  tried  as  initial  conditions  for  the 
plume.  In  addition,  a  constant  pressure  mixing  (fully  pressure 
equilibrated  to  ambient  pressure)  case  was  compared  to  the  more  precise 
overlaid  calculation  (mixing  layer  superposed  on  the  underexpanded 
plume  boundary).  In  this  particular  case,  there  were  no  significant 
differences  between  the  overlaid  and  constant  pressure  mixing  solutions 
because  the  degree  of  underexpansion  was  not  severe  (ratio  of  exit  to 
ambient  pressure  of  about  two). 

There  was  a  large  difference  between  the  predictions  using  the 
equilibrium  and  chemical  kinetic  initial  conditions.  The  kinetics 
solution  provided  a  higher  electron  number  density,  and  consequently 
the  attenuation  through  the  initial  region  of  the  plume  (along  those 
paths  penetrating  the  inviscid  core  of  the  plume)  was  too  large.  The 
difference  between  these  two  sets  of  initial  conditions  diminished 
downstream  as  the  chemistry  in  the  plume  mixing  region  began  to 
dominate  over  the  initial  concentrations.  However,  even  with  the  lower 
electron  concentrat Ions  of  the  equilibrium  initial  condition*,  the 
predicted  attenuation  was  too  large,  for  the  purpose  of  testing  PRF1C, 
we  took  the  approach  of  reducing  the  electron  concentration  uniformly 
throughout  the  plume  until  the  prediction*  and  measurements  matched  at 
the  first  measurement  location,  (about  z/rj  •  18)  near  the  end  of  the 
Inviscid  core.  By  reducing  the  electron  number  density  uniformly  bv  a 
multiplier  of  l/B,  It  was  possible  to  match  the  tnitlal  transverse 
attenuation  at  the  rhree  measurement  frequencies  of  9,  16,  and  12  GHz. 
These  results  ire  given  In  Figures  7,  8,  and  9,  for  these  three 
frequencies.  It  Is  Important  to  note  that  the  downstream  distribution 
of  attenuation  Is  also  closely  predicted.  By  calibrating  the  flowfteld 
at  the  first  measuremen'  station,  we  have  made  good  predictions  of  the 
attenuation  downstream  ,-s  well.  We  Infer  from  this  comparison  that, 
when  given  the  «*•  jper  flowfteld  description,  the  attenuation 
predictions  with  PRFIC  compare  favorably  with  the  three  sets  of 
measurements  at  different  frequencies. 

It  is  important  to  not*  that  the  plume  prediction  giving  the  best 
agreement  with  the  attenuation  measurements  wss  made  with  the 
two-equation  turbui»nc*  model,  without  the  compressibility  correction. 
The  flowfleld  with  the  compressibility  correction  gave  too  much 
attenuation  downstream.  The  exit  Nach  number  was  3.2,  and  the  Inviscid 
core  lengths  were  29rj  snd  l*rj  with  and  without  the  correction, 
respectively.  The  same  conclusion  was  rrsched  independent ly  in  s  study 
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FIGURE  3.  Transverse  Attenuation  at  16  GHz.  NWC  TP  5319,  Part  1, 
Case  22. 
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of  the  gas  dynamic  structure  of  plumes.10  This  sensitivity  suggests 
that  both  the  magnitude  and  distribution  of  plume  RF  interference 
properties  are  strongly  affected  by  our  ability  to  predict  turbulent 
flows . 

An  additional  set  of  transverse  attenuation  and  phase  predictions 
is  given  for  the  exhaust  plume  shown  in  Figures  3,  4,  and  5.  A  profile 
of  transverse  attenuation  from  the  axis  to  the  edge  of  the  plume  is 

given  in  Figure  10.  At  5  meters  downstream  (z/rj  ■  37),  the  peak 

electron  number  density  and  attenuation  are  located  in  an  annular 
volume  about  :ie  axis.  The  attenuation  shows  this,  as  the  peak  occurs 
at  a  radius  >  f  about  0.125  meters  ( ( y/ r j )  ~  1).  This  prediction 

suggests  that  the  peak  attenuation  does  not  necessarily  occur  for  a 
line-of  -sight  hrou^h  the  plume  axis,  at  least  for  those  regions  where 
afterburning  is  most  intense. 

Attenuation  predictions  using  the  ray  trace  option  are  also  shown 
in  Figure  10.  For  this  particular  example,  the  refracted  paths  yield 
attenuation  which  is  not  substantially  different  from  the  straight 

1  lne-of-sight .  When  there  is  a  difference,  the  attenuation  is 
uniformly  larger  for  the  refracted  path. 

The  phase  change  for  the  same  conditions  as  the  attenuation  in 

Figure  10  is  shown  in  Figure  11.  The  phase  difference,  , 


$ds  -  2i»s/X 


/X 


(44) 


between  that  along  the  actual  path  and  that  for  the  free  space 
1 ine-of-sight  is  shown.  The  peak  phase  shift  occurs  at  the  location  of 
peak  attenuation.  We  n.-te  that  the  coherent  RF  power  ratio  due  to  this 
phase  difference  is 


~*2iA* 


(45) 


The  phase  difference  predicted  for  thp  refracted  path  is  larger  by 
a  small  amount  in  the  highly  attenuating  region  for  this  case. 


TTTpe  rgaraent  ,  H.S.  ,  "Assessment  and  Recommendation  of  TWo-Equatlon 
Turbulence  Models  for  Rocket  and  Aircraft  Flume  Flowfleld 
Predictions,"  Naval  Weapons  Center,  China  Lake,  Calif.,  July  1982. 
(NWC  TP  6)64,  publication  UNCLASSIFIED). 
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DIAGONAL  ATTENUATION 

An  example  of  diagonal  attenuation  from  near  tail-on  observation 
angles  for  the  plumes  shown  in  Figures  3-5  is  given  in  Figure  12. 
These  predictions  show  the  rapid  onset  of  severe  attenuation  as  the 
line-of-sight  passes  through  the  edge  of  the  plume  and  into  the  region 
where  it  passes  through  nearly  the  entire  plume  length.  Calculations 
like  this  have  shown  that  the  attenuation  for  observation  on  the  same 
side  as  the  antenna  (g  >  0)  is  very  sensitive  to  the  radial  location  of 
the  antenna.  Additionally,  this  situation  is  one  in  which  diffraction 
by  the  plume  is  Important .  The  increased  attenuation  wi^h  decreasing 
angle  is  actually  less  than  that  shown  in  Figure  12,  when  the 
diffracted  wave  is  included.  For  this  particular  geometry,  the 
attenuation  due  to  absorption  calculated  by  PRFIC  is  a  worst  case,  and 
the  effects  of  diffraction  should  be  included  to  provide  an  accurate 
description  of  the  total  plume  interference. 

An  example  of  the  three-dimensional  capability  of  the  PRFIC  code  is 
given  in  Figure  13  in  which  the  attenuation  of  the  transmission  from  a 
rolling  missile  is  given.  These  calculations  are  for  a  large  Minuteman 
class  plume  at  5  km  altitude.  The  observation  is  from  an  aspect  angle 
of  150  degrees  (30  degrees  from  tail-on).  A  zero  roll  angle  occurs 
when  the  transmitting  and  receiving  antennas  are  on  opposite  sides  of 
the  plume.  The  attenuation  diminishes  as  the  antenna  rolls  out  of  the 
shadow  of  the  plume.  The  calculation  requires  the  capability  of  the 
code  to  deal  with  a  fully  three-dimensional  line-of-sight  through  the 
axisymmetric  plume.  This  observational  condition  is  also  one  in  which 
diffraction  effects  are  important. 


SCATTERING 

There  are  several  applications  of  the  calculation  of  scattering 
from  the  turbulent  fluctuations  in  the  plume  electrical  properties. 
Scattering  in  the  direction  of  propagation,  like  diffraction  by  the 
mean  flow,  increases  the  power  received  and  reduces  the  apparent 
attenuation.  In  the  limiting  case  where  there  is  no  direct  propagation 
path  from  transmitter  to  receiver,  scattering  and  diffraction  are  the 
only  means  of  signal  transmission.  In  the  case  of  pure  backscatter, 
scattering  by  turbulence  is  a  component  of  the  plume  radar  cross 
section  in  addition  to  that  scattered  by  the  mean  flow.  Finally,  the 
fluctuating  nature  of  scattering  by  the  turbulence,  and  the  fact  that 
it  has  a  large  Doppler  shifted  frequency  spectrum  because  of  the 
turbulent  velo  ity  distribution,  create  a  source  of  noise  that  is 
present  in  any  portion  of  the  received  signal  that  has  interacted  with 
the  exhaust  plume. 
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FIGURE  12.  Diagonal  Attenuation  at  Near  Tail-on  Aspect. 
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FIGURE  13.  Diagonal  Attenuation  at  Near  Tail 
Aspect  for  a  Rolling  Missile,  Minuteman  Class 
Plume  at  5  km  Altitude. 
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As  an  Initial  example  of  the  calculation  of  scattering,  we  include 
a  tabulation  of  the  scattered  power  for  the  transverse  attenuation 
conditions  shown  in  Figures  7,  8,  and  9  (see  Table  2).  This  is  the 
power  added  back  into  the  direction  of  propagation  that  reduces  the 
effective  attenuation.  The  contribution  is  that  added  incoherently, 
without  regard  to  phase,  and  does  not  include  attenuation  of  the 
incident  and  scattered  wave  (the  added  effect  of  attenuation  is  small 
for  this  plume.)  In  all  cases,  the  power  received  by  scattering  is 
very  much  less  chan  that  absorbed  from  the  directly  transmitted  wave. 


TABLE  2.  Transverse  Attenuation  and  Scattering 

Predictions  for  cases  21,  12,  23,  of  NWC  TP  5319,  Part  1  (Reference  1) 
(scattering  calculations  do  not  include  attenuation) 


9  GHz  1 

16 

GHz 

32 

GHz 

Attn.dB  | 

Scatt  ,dB 

Attn.dB 

'  Scatt, dB 

17.7 

1.34 

!  -91.79 

1.104 

-94.34 

0.553 

-100.16 

35.4 

t 

2.69 

-87.06 

1.62  1 

1 

-88.85 

0.599 

-  94.85 

53. 1 

2. 19 

-86.21 

1.31 

-89.28 

0.470 

-  97.27 

70.7 

1.03 

-89.07 

0.608 

-92.80 

{ 

0.  169 

-100.36 

88.4 

0.373 

-95.74 

i  1 

0.203 

-99.98 

0.015 

-107.37 

These  scattering  results  were  computed  with  uniform  gain  for  both 
the  transmitter  and  receiver.  If  a  highly  directional  antenna  is  used 
(for  example,  with  a  gain  of  20  dB)  for  both  the  transmitter  and 
receiver,  the  scattering  contribution  is  still  small  compared  to  that 
of  the  direct  transmission  for  this  particular  plume. 

The  contribution  to  the  received  signal  due  only  to  scattering  was 
computed  for  a  fixed  receiver  location  and  direction  (a  radial  location 
of  y/ r j  -  1.86,  pointed  downstream,  parallel  to  the  plume  axis).  The 
plume  properties  are  those  shown  in  Figures  3-5.  The  scattered  power 
received  as  the  transmitter  was  moved  around  the  vehicle  is  shown  in 
Figure  14.  Cases  both  with  and  without  attenuation  of  the  scattered 
waves  are  shown.  The  polar  plot  is  for  the  situation  in  which  the 
receiver,  transmitter,  and  plume  axis  are  all  in  the  same  plane.  This 
result  shows  that,  for  this  particular  plume,  the  effect  of  attenua¬ 
tion  of  the  scattered  power  contributes  about  a  10  dB  loss  to  that 
scattered  without  attenuation.  In  both  cases,  there  is  ver;  little 
aspect  angle  dependence  in  the  forward  hemisphere  of  scattering 


/ 
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(0  <  a  <  90  degrees  and  270  <  a  <  360  degrees).  There  is  an  aspect 
angle  dependence  In  the  rearward  hemisphere,  and  there  is  an  asymmetry 

about  the  plume  axis  because  the  transmitter  and  receiver  are  on  the 

same  side  of  the  plume  for  180  degrees  <  a  <  360  degrees  and  on 

opposite  sides  for  0  *  a  <  180  degrees.  It  is  Important  to  note  that, 
when  the  aspect  angle  is  in  the  forward  hemisphere,  there  is  no  direct 
transmission,  and  the  received  power  is  only  that  scattered  by  the  mean 
flow  and  the  turbulence. 

Sample  calculations  of  the  aspect  dependent  backscatter  (i.jdar 
cross  section)  by  the  turbulent  plume  flow  are  shown  in  Figure  15.  The 
same  plume  flowfield  was  used,  and  results  for  scattering  with  and 

without  attenuation  are  shown.  In  the  forward  hemisphere,  no  allowance 
was  made  for  occultatlon  by  the  missile  body.  These  results  are 
symmetric  with  respect  to  the  plume  axis.  There  is  very  little  aspect 
angle  dependence  of  tne  scattering  without  attenuation  (at  sufficiently 
large  distances,  the  total  plume  cross  section  depends  only  on  the 
number  of  electrons  and  their  fluctuations).  In  contrast,  the 
scattering  with  attenuation  decreases  sharply  at  the  ne;  r  tai L-on 
aspect  angle  because  of  absorption  within  the  larger  path  lengths 
through  the  plume. 


DOPPLER  SHIFT 

The  Doppler  shift  of  the  scattered  power  due  both  to  the  mean 
relative  missile,  plume,  and  receiver  velocities  and  to  the  turbulent 
fluctuations  contributes  a  frequency  spectrum  of  power.  An  example  of 
this  spectrum  is  given  in  Figure  16  for  a  Minuteman  class  missile  with 
a  velocity  of  503  m/s.  These  calculations  are  for  the  backseat tered 
power,  and  the  spectra  are  shewn  with  and  without  the  effect  of 
attenuation  of  Che  scattered  power.  There  is  an  important  effect  of 
the  attenuation  on  the  frequency  spectrum;  and,  qualitatively  at 
least,  the  spectra  with  attenuation  have  the  general  appearance  ot  the 
measurement  of  AM  noise  Reference  2). 

Additional  examples  of  the  scattering  frequency  spectra  are  shown 
in  Figures  17  and  18,  which  are  for  aspect  angles  of  45  and 
135  degrees,  respectively.  These  two  observations  correspond  to  the 
missile  approaching  the  observer  and  receding  from  the  observer.  The 
two  results  show  the  change  in  spectra  for  these  two  relative 
velocities.  The  static  case  is  superposed  for  comparison. 

Finally,  we  show  several  frequency  distributions  of  the  scattered 
power  for  different  aspect  angles  in  Figure  19.  These  results  are  for 
the  Minuteman  class  plume,  and  are  typical  AM  noise  spectra.  Although 
no  quantitative  comparisons  have  yet  been  made,  it  is  encouraging  that 
these  predicted  spectra  have  the  general  appearance  of  the  several 
measurements  summarized  in  Reference  2. 
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FIGURE  15.  Plume  Radar  Cross  Section  (plume  properties  of 
Figures  3-5) . 


37 


SM2-2#  dVBV 


NWC  TP  6386 


Aspect  Angle  =  135  Deg 
*0  « 9  GHz 


'o  ,  *Hz 

FIGURE  16.  Ooppler  Shifted  Frequency  Spectra  for  Backscatter 
from  a  Mnutoman  Claes  Plume  at  1  .5  km  Altitude,  Velocity  of 
503  m/i ,  135  deqrees  Aspect  Anqle. 
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tO 
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FIGURE  18.  Doppler  Shifted  Frequency  Spectra  for 
Backscatter  from  a  fllnuteman  Class  Plume  at  1.5  km 
Altitude,  Velocity  of  503  m/s.  Aspect  Angle  of 
1 35  degrees . 
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FIGURE  19.  Frequency  Spectra  for  Backscattering  from  a  Minuteman 
Class  Plume  at  1.5  km  Altitude,  Velocity  of  503  m/s. 
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These  are  representative  samples  of  the  Doppler  spectra  of  the 
scattered  signal.  Their  shape  is  affected  to  some  extent  by  the  choice 
of  a  Gaussian  probability  density  for  the  magnitude  of  the  turbulent 
velocity  fluctuations.  Other  forms  have  not  been  examined,  and  so  the 
sensitivity  of  the  spectrum  to  this  choice  is  unknown.  However,  it  is 
an  easy  task  to  change  the  probability  density  in  the  code,  and  an 
assessment  of  various  choices  against  experimental  data  can  easily  be 
made.  We  anticipate  that  PRFIC  will  be  a  valuable  tool  for  this  kind 
of  analysis  and  predictions. 


SUMMARY  OF  CAPABILITIES  AND  LIMITATIONS 


The  primary  capabilities  and  limitations  of  PRFIC  are  summarized 
below. 

1.  Flowfield  Input  -  SPF  format:  radial  distribution  of  properties  at 
successive  axial  stations 

2.  Geometry  -  Arbitrary  location,  orientation,  and  velocity  of  the 
missile,  receiver,  and  transmitter  (see  Figure  6) 

(This  generality  was  included  to  account  for  transmission  between  a 
moving  missile  and  a  fixed  or  moving  observer,  static  test 
situations,  and  radar  cross  section  predictions  with  a  single 
code. ) 

3.  Direct  Transmission 

a.  Attenuation  and  phase  predictions,  Doppler  shift 

b.  Llne-of-sight ,  or  ray  trace  option 

(The  ray  trace  calculation  is  not  reliable  in  its  current  form 
if  there  are  very  strong  gradients  in  index  of  refraction.) 

4.  Scattering 

a.  Single  scattering  from  turbulence-induced  index  of  refraction 
fluctuations,  with  attenuation  of  the  Incident  and  scattered 
wave 

b.  Doppler  frequency  spectra  of  the  scattered  power 

(The  current  code  uses  a  Gaussian  distribution  of  the 

fluctuating  velocity  probability  density.) 

c.  Electron  number  density  fluctuations  which  are  obtained  from 
the  SPF 

(The  modeling  of  n^  is  confined  to  the  SPF,  the  specification 
of  ,  given  n^,  is  in  PRFIC.) 

(The  current  version  of  PRFIC  keeps  track  of  the  amplitude  of 
the  frequency  spectra,  the  phase  of  the  scattered  wave  is 
computed,  but  not  currently  utilized.) 

(The  scattering  path  can  be  computed  with  either  a 
1  ine-ol -s  i  ght  <>r  r.iv  trace  option.) 
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SUGGESTED  IMPROVEMENTS 


There  are  two  classes  of  improvements  to  PRFIC  that  can  be  easily 
made.  One  set  of  improvements  involves  the  addition  of  new  physical 
phenomena  to  the  code.  In  this  category,  the  authors  suggest  the 
following : 

1.  Focusing  effects  of  the  refraction 

2.  Surface  scattering  effects  (diffraction)  at  the  overdense 
surface 

Simple  diffraction  models  have  been  quite  successful  in  predicting 
attenuation  by  plumes  with  moderate  to  high  electron  density  levels 
(References  1  and  3).  However  such  models  do  not  contain  the 
scientific  principles  or  details  needed  to  assure  accurate  predictions 
in  all  cases.  As  originally  conceived,  PRFIC  was  expected  to  give  the 
same  effects  as  a  plume  diffraction  calculation  if  absorption, 
scattering,  and  phase  shift  calculations  are  combined  and  a 
sufficiently  fine  computational  grid  is  used.  If  PRFIC,  as  currently 
implemented,  has  a  major  weakness,  it  is  probably  the  lack  of  a 
mechanism  to  account  for  surface  scattering  at  the  boundary  that 
defines  the  overdense  surface  of  plumes  with  high  electron  densities. 
The  results  of  scattering  in  that  region  should  be  quite  similar  to 
those  of  a  diffraction  calculation. 

The  authors  suggest  that  PRFIC  is  appropriate  to  serve  as  a  basis 
for  a  solution  which  Incorporates  diffraction  by  the  mean  flcwfield 
and/or  scattering  by  turbulence  at  the  boundary  of  the  overdense 
region,  whichever  is  determined  to  be  the  more  appropriate  approach. 
The  addition  of  a  diffraction  capability  has  the  added  advantage  that 
it  would  provide  the  code  with  the  ability  to  calculate  missile  body 
effects  on  the  radiation  field  in  cases  where  body  mounted  antennas  are 
involved.  These  extensions  to  PRFIC,  along  with  the  current 
capabilities,  would  provide  an  entire  plume  interference  computational 
capability  in  one  unified  code. 

Another  category  of  suggested  improvements  involves  strictly  the 
computational  aspects  of  the  current  code.  Some  improvement  in  the 
numerical  solution  for  the  refracted  ray  path  is  needed  to  make  it  more 
reliable.  If  there  are  very  severe  gradients  in  n,  the  current 
integration  scheme  does  not  always  converge.  In  general,  the  solution 
converges  for  tr. Rectories  that  remain  in  the  same  plane.  For  example, 
convergence  is  obtained  for  ray  paths  that  remain  either  nearly 
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perpendicular  to  the  nozzle  xis  (transverse  attenuation),  or  nearly 
parallel  to  the  axis.  Solutions  have  been  obtained  for  these  cases 
even  for  the  most  severe  index  of  refraction  gradients.  The  solution 
seems  to  be  most  sensitive  along  those  trajectories  which  do  not  remain 
planar.  Our  Initial  attempts  to  define  the  criteria  for  convergent 
solutions  more  precisely  have  not  been  successful.  Since  the  departure 
of  the  ray  trajectory  from  the  straight  line-of-sight  depends  both  upon 
the  index  of  refraction  gradient  normal  to  the  trajectory  and  on  the 
path  length,  the  probability  of  non-convergence  increases  with  the 
integral  of  the  gradient  along  the  ray  path.  However,  at  this  point, 
we  do  not  have  a  specific  criterion  to  define  those  conditions  where 
convergence  is  assured. 

The  two-point  boundary  value  problem  is  one  of  the  more  difficult 
numerical  integrations  to  perform,  especially  when  the  solution  is  very 
sensitive  to  the  initial  conditions  as  is  the  case  in  the  present 
application.  The  algorithm  used  in  the  current  version  of  PRFIC  can  be 
improved  with  the  addition  of  a  more  refined  logic  to  improve 
convergence  for  the  more  difficult  ray  trajectories. 

Finally,  an  improvement  in  the  integration  procedure  for  both  Che 
straight  line-of-sight  and  refraction  solution  is  suggested.  The 
current  procedure  is  to  choose  the  integration  step  size  ja  priori  and 
to  use  uniform  increments  throughout  the  entire  path  length  within  the 
plume  velocity  boundaries.  However,  often  the  significant  electron 
number  density  is  confined  to  a  smaller,  electrical  plume.  The  small 
integration  steps  outside  the  effective  electrical  plume  are  wasted. 
An  improvement  would  be  to  let  the  step  size  vary  throughout  the  plume, 
and  depend  on  the  gradient  of  the  plume  electrical  properties  along  the 
line-of-sight.  This  is  an  easily  implemented  improvement  that  would 
yield  accurate  solutions  with  less  computational  time. 


-UMMARY  AND  CONCLUSIONS 


We  have  presented  a  summary  of  the  physics  of  the  interaction  of 
electromagnetic  waves  with  radar  frequencies  with  a  weakly  ionized 
rocket  exhaust  plume.  A  computer  code,  the  Naval  Weapons  Center  Flume 
Radar  Frequency  Interference  Code  (PRFIC),  has  been  developed  to 
calculate  several  aspects  of  this  interaction.  The  code  relies  on  the 
description  of  the  exhaust  plume  givpn  by  the  JANNAF  Standard  Plume 
Flowfield  Model  (SPF).  Its  current  capabilities  in  addition  to  the  SPF 
plume  structure  are: 

1.  full  three-dimensional  geometry  of  the  transmitter,  mlssil", 
and  receiver  locations 

attenuation  predictions  for  the  directly  transmitted  power 
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along  line-of-slght  or  refracted  ray  trajectories 

3.  scattering  from  turbulent  fluctuations  of  electron  number 
density 

4.  Doppler  shift  of  the  directly  transmitted  and  scattered 
signal 

There  are  two  additional  aspects  of  the  refracted  ray  that  are  not 
currently  included  in  PRFIC.  One  is  the  multipath  effect  in  which  rays 
leaving  the  transmitter  in  different  directions  can  arrive  at  the 
receiver.  The  multiple  effect  can  be  obtained  from  PRFIC  if  enough 
raypaths  are  calculated.  However,  it  is  not  an  automated  part  of  the 
calculation  sequence.  The  second  effect  is  that  of  focusing  the 
incident  energy  (see  Appendix  C).  Both  of  these  effects  are  directly 
related  to  the  refractive  properties  of  the  plume  and  are  absent  in  a 
1 ine-of-sight  approximation.  They  are  not  included  in  the  present, 
first  version  of  the  code.  We  have  no  a  priori  indication  of  their 
importance  in  actual  exhaust  plume  attenuation  predictions.  Subsequent 
versions  of  PRFIC  should  account  for  these  additional  effects  and 
calculations  should  be  made  to  determine  the  circumstances  for  which 
they  are  important. 

This  code  has  the  potential  for  dealing  with  many  of  the  variables 
that  affect  plume  RF  interference.  With  the  SPF,  we  can  assess  the 
effects  of  propellant  composition,  vehicle  size  (thrust),  and  flight 
conditions.  With  PRFIC,  we  can  predict  the  effect  of  the  geometry  of 
transmission  (antenna  locations),  frequency,  and  turbulence  properties. 
Combined,  these  capabilities  allow  a  very  complete  description  of 
several  aspects  of  plume  RF  interference. 

We  have  given  some  examples  of  the  predictions  by  both  the  SPF  and 
PRFIC.  These  are  included  to  illustrate  their  capabilities  and 
limitations.  Improvements  have  been  identified.  However,  the  full 
capability  remains  to  be  verified  by  intensive  comparisons  with 
detailed  measurements,  which  we  strongly  recommend. 
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Appendix  A 
USER'S  MANUAL 


INTRODUCTION 

This  appendix  is  the  user's  manual  for  the  Naval  Weapons  Center 
Plume  Radar  Frequency  Interference  Code  (PRFIC).  It  documents  the  code 
structure  and  detailed  input  instructions.  Recommended  values  for 
those  variables  that  define  convergence  criteria  and  number  of 
iterations  are  given.  Sample  input  and  results  from  test  cases  are 
included.  The  physics  involved  in  the  calculations  are  described  in 
the  main  text  of  this  report,  and  the  detailed  software  description  is 
given  in  Appendix  B. 


CODE  STRUCTURE 

Essentially  each  type  of  calculation  performed  in  PRFIC  is 
separated  into  individual  subroutines.  A  flow  chart  that  identifies 
the  subroutines,  and  the  sequence  in  which  they  are  called,  is  given  in 
Figure  A-l.  A  complete  list  of  subroutines  is  given  in  Table  A-l. 
More  complete  descriptions  of  the  most  important  subroutines  are  given 
in  Appendix  B.  The  input  flowfield  file  with  the  appropriate  variables 
defining  the  electrical  properties  of  the  plume  is  obtained  from  the 
JANNAF  SPF  code  or  an  equivalent  calculation. 


Table  A-l.  List  of  Subroutines  and  Their  Function  in  PRFIC 


Subroutine 

Function 

PRFIAG: 

Computes 

antenna  gain 

PRFICM: 

Common 

PRFIDT : 

Sets  up  initial  and  final  points  for  integration  through 
plume 

PRFIEP: 

Computes  electrical  properties 

PRFIFF : 

Reads  flowfield 

PRFIIC : 

Reads  changes  for  multiple  runs 

PRFIIN: 

Reads  input 

PRFIIT: 

Performs  integration  through  plume 
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Table  A-l.  (Continued) 

Subroutine  Function 

PRFILI:  Linear  interpolation 

PRFIMM:  Main  calling  subroutine 

PRFIPC:  Computes  attenuation  and  phase  constants 

PRFIPl:  Computes  intersection  of  line-of-sight  with  plume 

PRFIPL:  Printer  plot 

PRFIRZ:  Computes  plume  radius 

PRFISC:  Performs  volume  Integration  for  scattering 

DUMP  :  Writes  input  conditions 

PDF  :  Turbulent  velocity  probability  density  function 

(currently  Gaussian) 

SCAT  :  Computes  total  scattering  cross  section 

SCALP  :  Computes  scalar  product  of  two  vectors 

SFVMV  :  Exchanges  one-dimensional  arrays 

(system  subroutine) 

CODE  INPUT 

A  description  of  each  input  variable  and  the  Instructions  for 
running  multiple  cases  are  given  in  this  section.  The  coordinate 
system  used  in  the  code  was  shown  in  Figure  6  of  this  report.  The 
primary  input  is  listed  in  Table  A-2. 
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/ 


FIGURE  A-l .  Structure  of  the  Naval  Weapons  Center  Flume 
Radar  frequency  interference  Code  (PRFIC). 
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Table  A-2.  Primary  Input 


Card 

Columns 

Fo  rma  t 

Var iabl 

1 

1-20 

20A1 

- 

2 

1-72 

18A4 

- 

3 

1-10* 

E10.3 

1 

11-20 

no. 4 

2 

21-30  F10.4  3 


4 


1-10 

F10.4 

1 

11-20 

F10.4 

2 

21-30 

F10.4 

3 

No .  Variable 

Flowfield  file  name, 
including  extension 
(i.e.,  FILENAME.DAT) 

Title 

BKREQ,  radar 
frequency,  Hz 

ZLL,  upstream  extent 
of  flowfield  to  be 
considered . 

ZLL  -  0,  the  first 
station  on  flowfield 
file  is  used. 

ZLL  >  0,  only  Z  >  ZLL 
considered . 

ZUL,  downstream  extent 
of  flowfield  to  be 
considered . 

ZUL  <  0,  the  last 
station  on  the  flowfield 
file  is  used. 

ZUL  >  0,  only 
Z  <  ZUL  considered. 

XT(  1 ) , XT(  2  )  ,  XT(  3 )  , 
components  of  position 
vector  to  the  radar 
transmitter,  m. 

XT(1)  -  x-component 

XT(2)  -  y-coraponent 

XT(3)  -  z-component 


*In  this  and  all  subsequent  cards,  a  list-directed  input  can  be  used  in 
which  the  format  is  ignored  and  the  variables  are  input  in  order  and 
separated  by  commas. 
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Table  A-2.  (Continued) 


Card  Columns  Format  Variable  No. 


Variable 


6 


1-10  F10.4  1 


1-10  F10.4  2 


1-10  F10.4  3 


6 


1-10  F10. 4  1 


DCT(l),  DCT( 2) ,  DCT( 3)  ,♦ 
Direction  of  radar 
transmitter  antenna 
with  respect  to  the 
x,y,z  coordinate 
system,  degrees 

DCT( 1 ) ,  angle  between 
antenna  direction  and 
x-axis 

DCT(2),  angle  between 
antenna  direction  and 
y-axis 

DCT(3),  angle  between 
antenna  direction  and 
t-axis 

DCT0(1),  DCT0( 2 )  ,** 
DCTO(3),  Directions  of 
a  reference  at  the 
transmitter  antenna, 
degrees 

DCT0( 1 ) ,  angle  between 
reference  direction  and 
x-axis 


♦See  Figure  ft  where  these  angles  are  called  a,  S,  and  t,  respectively. 
Subscripts,  T,  R,  and  M,  refer  to  the  transmitter,  receiver,  and 
missile.  For  all  sets  of  these  directions,  it  is  necessary  that 
coa^a  +  cos^S  +  cos^y  ■  1, 

♦♦These  reference  directions,  for  the  transmitter,  receiver,  and 
missile,  are  included  for  future  use  (polarization  effects,  for 
example);  none  are  used  explicitly  in  the  present  calculations. 
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Table  A-2.  (Continued) 


Card 

Co  lunllis 

Format 

Variable  No. 

Variable 

11-10 

F10.4 

2 

DCT0{2),  angle  between 
reference  direction  and 
y-axls 

21-20 

F10.4 

1 

DCT0(3),  angle  between 
reference  direction  and 
z-axls 

7 

vr(l),  VT(2),  VT( 3) , 
components  of  velocity 
vector  of  the  radar 
transmitter,  m/s 

1-10 

F10.4 

1 

VT(1)  component  in 
x-directlon 

11-20 

F10.4 

•> 

VT(2)  component  in 
y-dlrectlon 

21-30 

FI0.4 

3 

VT(3)  component  in 
z-dlrectlon 

S 

XR(  1 ) ,  XR( 2 )  •  XRO) 
Component  of  position 
vector  to  the  receiver, 

1-10 

FI  0.4 

1 

XR( 1 ) ,  x-conponent 

11-20 

F10.4 

2 

XR(2),  y-co«ponent 

21-30 

F10.4 

3 

XR(1),  *-component 
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Table  A-2.  (Continued) 


Card 

Columns 

Format 

Variable  No. 

Variable 

9 

DCR(l),  DCR( 2 ) ,  DCR( 3) , 
Directions  of  the  radar 
receiver  antenna,  degrees 

1-10 

F10.4 

1 

DCR(l),  angle  between 
receiver  direction 
and  x-axis 

11-20 

F10.4 

2 

DCR(2),  angle  between 
receiver  direction 
and  y-axis 

21-30 

F10.4 

3 

DCR(3),  angle  between 
receiver  direction 
and  Z-axis 

in 

DCRO(l),  DCR0( 2 ) ,  DCR0(3) 

1-10 

F10. 4 

1 

DCT0( 1 ) ,  angle  between 
reference  direction 
and  x-axis 

11-20 

F10. 4 

2 

DCT0(2),  angle  between 
reference  direction 
and  y-axis 

21-30 

F10. 4 

3 

DCT0(3),  angle  between 
reference  direction 
and  z-axls 

11 

VR(1),  VR(2),  VR(  3) 
Components  of  velocity 
vector  of  radar 
receiver,  a/s 

1-10 

F10. 4 

1 

VR( 1 ) ,  x-componcnt 

11-20 

F10. 4 

2 

VR(2),  y-coraponent 

21-10 

F10.4 

3 

VR(3),  z-coraponent 
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Table  A-2.  (Continued) 


Card 

Col umns 

Fo  rma  t 

Variable  No. 

Variable 

12 

XM(  1 ) ,  XM(  2 ) ,  XK(  3) 

Components  of  position 
vector  of  missile 
exit  plane,  ra 

1-10 

F10.4 

1 

XM(1),  x-coaponent 

11-20 

F10.4 

2 

XM(2),  y-component 

21-10 

F10.4 

3 

XM(3),  z-coraponent 

Note:  Normally,  the 
missile  exit  plane 
is  at  the  origin 
of  the  coordinate 
system  so  that 

XM(  1 )  -  XM(  2)  -  XM(  3)  -  0. 

13 

DCM(l),  DCM( 2 ) ,  DCM( 3) 

Angles  of  normal  to  the 
missile  exit  plane, 
degrees 

1-10 

F10.4 

1 

DCM(l),  angle  between 
normal  and  x-axis 

1 1-20 

F10.4 

2 

DCM(2),  angle  between 
normal  and  y-axis 

21-30 

Fl(). 4 

3 

DCM(3),  angle  between 
normal  and  z-axls 

Note:  Normally,  the 
missile  axis  Is 
pointed  along  the 
z-axls,  so  that 

DCM(l)  -  DCM( 2 )  -  90  degrees, 
DCM(3)  -  0. 

14 

DCMO(I),  DCMO( 2 ) ,  DCM0( 3) 
Angles  of  a  reference 
direction  at  the  missile 
exit  plume,  degrees 
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Table  A-2.  (Continued) 


Columns 

Format 

Variable  No. 

Variable 

1-10 

FJ0.4 

1 

DCMO(l),  angle  between 
reference  and  x-axis 

11-20 

F10.4 

2 

DCM0(2),  angle  between 
reference  and  y-axis 

21-30 

F10.4 

3 

DCM0(3),  angle  between 
reference  and  z-axls 

VM( 1 )  ,  VM( 2)  ,  VM( 3)  , 
components  of  missile 
velocity  vector,  m/s 

1-10 

F10.4 

1 

VM( 1 ) ,  component  in 
x-direction 

11-10 

F10.4 

2 

VM(2),  component  in 
y-direction 

21-30 

F10.4 

2 

VM( 3)  ,  component  in 
z-direction 

1-10 

F10.4 

1 

NCUTS  -  number  of 
steps  used  in 
integration  for 
attenuation  along 
the  ray  path  through 
the  plume 
( NCUTS  -  10-40) 

11-20 

F10.4 

2 

RC  -  overrelaxat Ion 
factor  in  ray  trace 
iteration 
(0.5  <  RC  <  1.0) 

21-30 

110.4 

3 

NIT  -  number  of 
iterations  allowed 
in  ray  trace 
calcul at  Ion 
(10  <  NIT  <  20) 
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Table  A-2.  (Continued) 


?ard  Columns  Format 
17 


1-10  F10.4 

U-20  F10. 4 


Variable  No.  Variable 

EPSA(l),  EPSA( 2) , 
error  criteria  for 
convergence  of 
ray  trace  Iterations,  m 

'  EPSA( 1 )  —  convergence 

in  the  interior  of 
the  plume 

2  EPSA(2)  -  convergence 

at  the  receiving 
antenna 


18 


EPSA(2)  <  radius  of 
antenna 


RAYTR  -  ray  trace 
option  flag. 

RAYTR  -  T,  ray 
trace  calculation 
RAYTR  ■  F,  line  of 
sight  calculation 

FSCAT  -  scattering 
calculation  flag. 
FSCAT  »  T,  radar 
cross  section 
calculation 
FSCAT  ■  F,  bypass 
radar  cross  section 
calculation 

Note:  If  RAYTR  -  X 
and  FSCAT  -  T,  the 
radar  cross  section 
calculation  Is  also 
done  with  the  ray 
trace  option. 
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Table  A-2.  (Continued) 


Card  Columns  Format  Variable  No.  Variable 

The  following  arrays  specifying  antenna  gains  are  not 
required  if  FSCAT  -  F. 

19 

1-5  15  1  NANG(l)  -  Number 

of  angles  defining 
transmitter  gain 
distribution. 

NANG(l)  >  0, 
read  NANG(l)  points 
NANG(l)  *  0,  uniform 
gain,  0-90  degreees 
GAIN( 1 , J)  -  0. 

NANG(l)  <  0, 
transmitter  gain 
set  equal  to 
receiver  gain. 

NANG(l)  <  50 

6-10  15  2  NANG( 2)  -  number  of 

angles  defining 
receiver  gain 
distribution 
NANG( 2)  >  0  read 
NANG( 2 )  points. 

NANG( 2)  -  0 
uniform  gain. 

0-90  degrees 
GAIN( 2 , J)  -  0. 

NANG( 2)  <  0 
receiver  gain 
set  equal  to 
transmitter  gain 
NANG(2)  <  50 

20  1-80  10F8.0  ANGLE(l.J),  J  -  1,  NANG( 1 ) , 

Angles  for 

transmitter  gain,  degrees 
in  ascending  order 

20+  1-80  10F8.0 
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Table  A-2.  (Continued) 


Card  Columns  Format  Variable  No.  Variable 

20+  1-80  10F8.0  ANGLEC  2 , J)  J  «  1,  NANG( 2) , 

Angles  for  receiver 

gain,  degrees  In  ascending 

order 

20+  1-80  10F8.0  CAIN( 2 , J) ,  J  -  1,  NANG( 2 ) , 

Gain  for  receiver,  dB 


MULTIPLE  CASES 

The  code  can  accommodate  multiple  cases  for  the  same  flowfiold  file 
and  antenna  gain  patterns.  That  is,  multiple  cases  for  which  changes 
occur  in  cards  3-18  (excluding  BFREQ.  which  occurs  on  card  3)  can  be 
run.  Only  those  variables  which  change  need  be  specified.  The  format 
is  given  below  and  also  in  the  FORTRAN  Source  Code  (in  PRFIIN). 

After  the  entire  set  of  input  for  the  first  card  (including  antenna 
gains) : 


Card 

1 

2 


2 


Columns  Format 

1  free  format 

1  free  format 

IC,  IV,  CV,... 


Variable 


NC  -  total  number  of 
changes  for  next  run 

IC  *  card  number  of 
change  3  *  IC  <  18 
IV  *  number  of  variables  on 
card  IC  to  be  changed 
(excluding  IV* l,  if  IC»3) 

CV  *  new  value  of  variable 

repeat  the  sequence  IC,  IV, 
CV  on  card  2  or  following 
cards  up  to  NC  times. 


3 
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This  sequence  is  repeated  for  each  subsequent  run.  However,  once 
the  changes  are  made,  they  remain  in  effect  until  changed  again.  Also 
note  that  if  an  antenna  gain  is  used,  it  must  be  read  in  the  initial 
data  set.  The  subroutine  that  reads  subsequent  cases  does  not  go  back 
to  read  the  gain  pattern. 


OUTPUT 

Page  1.  Date,  title,  and  all  input  quantities.  Geometric  and  velocity 
inputs  are  grouped  for  the  transmitter,  receiver,  and  missile. 

Page  2.  Local  and  cumulative  properties  for  beam  propagation  through 
the  plume.  If  the  ray  trace  option  is  used,  only  the  final, 
converged  path  is  printed.  The  columns  are: 


X.Y.Z 


ALPHA,  DE 

BETA*L 

N 


Position  of  points  along  the  ray  path,  in  meters 
beginning  at  the  transmitter  and  ending  at  the 
receiver. 

Cumulative  attenuation  of  beam  power,  in  dB. 


-Jf 


:  Cumulative  phase  along  path: 

:  Index  of  refraction. 

ALPHA,  1/M:  Local  attenuation  constant,  m-1. 

BETA,  1/M  :  Local  phase  constant,  m-1. 

K*uAMBDA  :  Local  product  of  wavenumber 
macroscale . 


0  ds  . 


and  turbulence 


SIGMA,  1/M:  Local  scattering  parameter,  m-1  (that  portion 
independent  of  angle). 


(0.0528)8” J  rTp7n2  n|  rI  A* 
T 


!  l  +  (ver)/o))2  ] 


U,  M/S  :  Local  axial  velocity,  m/s. 

(Q/U)**2  :  Local  turbulence  intensity,  normalized  by  local 

axial  velocity. 

These  quantities  are  printed  for  the  beginning  and  end  of  the 
ray  path  and  for  the  NCUTS-1  points  within  the  plume. 
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Note:  If  there  is  no  intersection  of  the  ray  path  with  the 
plume,  or  if  the  angle  between  the  ray  and  normal  to 
either  the  transmitter  or  receiver  antenna  is  greater 
than  90  degrees,  a  flag  will  be  printed  and  no 
Integration  will  be  performed. 

A  summary  of  the  cumulative  attenuation,  in  dB,  and 
phase  for  the  entire  path  is  given  (even  if  there  is  no 
Intersection  with  the  plume).  Also,  a  phase  difference 
(actual  phase  -  2irs/.\),  where  s  is  the  geometric  length 
between  the  transmitter  and  receiver,  is  given.  The 
Doppler  shifted  frequency  of  the  transmitted  ray  is 
printed,  in  kHz. 

Page  3.  Listing  of  transmitter  and  receiver  gain  patterns  (same  as 
input)  if  the  scattering  option  is  chosen. 

Page  4.  Frequency  spectrum  of  scattered  power.  The  columns  are  the 
number  of  points;  the  difference  between  the  shifted 
frequency  and  transmitter  frequency,  in  kHz  ((v  -  vo)/i000); 
a  normalized  frequency  shift,  (v  -  v0)c/v0,  and  the  scattered 
power,  in  dB/Hz,  referenced  to  the  total  transmitter  power. 

Frequency  integrated  scattered  power,  dB.  The  integral  with 
respect  to  frequency  of  the  spectral  power. 

Page  5.  A  printer  plot  of  the  frequency  spectrum,  dB/Hz,  vs.  v  -  vQ , 
in  kHz. 


DIAGNOSTICS 

There  are  several  built-in  diagnostic  messages  which  are  printed. 
They  are  identified  below,  along  with  the  subroutine  in  which  they  are 
located. 

Subroutine  Message/Meaning/Program  Action 

PXFIFF :  FLOWFIELD  E'JEEDS  RADIAL  DIMENSION 

There  are  ’.ore  than  50  radial  points  in  the  flowfield  file; 
program  °’.ops. 

PRFIDT :  CO-LOCATED  ANTENNAS 

There  is  no  attenuation  path  between  the  transmitter  and 
receiver;  computation  continues  if  FSCAT  »  T. 


* 
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NO  DIRECT  PATH  TO  RECEIVER,  S*NT  <  0 

The  ray  path  is  at  an  angle  of  more  than  90  degrees  with 

respect  to  the  normal  to  the  transmitter  antenna;  computation 
continues  if  FSCAT  «  T. 

NO  direct  path  to  receiver,  s*nr  <  0 

The  ray  path  is  at  an  angle  of  more  than  90  degrees  with 

respect  to  the  normal  to  the  receiver  antenna;  computation 
continues  if  FSCAT  =*  T. 

IT  -  MAX,  NO  CONVERGENCE 

The  ray  trace  calculation  failed  to  converge  to  an 
error  <  EPSA*(1)  in  NIT  iterations;  program  stops. 

PRFIPC :  PLUME  IS  LOCALLY  OVERDENSE 

(w p/oj ) / ( 1  +  (ven/ji)2)  at  some  point;  computation  continues. 

COMMON  PROBLEMS 

\ 

It  is  anticipated  that  most  problems  will  arise  in  the  ray  trace 
option.  If  a  portion  of  the  plume  is  overdense,  there  are  no  rays  that 

can  propagate  through  that  region,  and  a  ray  trace  calculation  will  not 

converge.  The  method  of  successively  computing  the  curved  path  from 
the  transmitter  towards  the  receiver  until  it  intersects  within  a 
prescribed  distance  EPSA(2)  proved  to  require  many  iterations  for  some 
paths.  The  solution  requires  experience  and  may  depend  on  relaxing  the 
error  criterion.  The  most  sensitive  paths  were  found  to  be  those 
requiring  corrections  in  two  directions  (i.e.,  those  paths  that  were 
not  wholly  contained  in  the  plane  of  symmetry  of  the  plume). 


LOGICAL  UNIT  ASSIGNMENTS 

The  following  logical  unit  (LUN)  assignments  are  used: 


LUN  Assignment 

5  Input 

6  Output 

7  Flowfield  file 
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SAMPLE  CASES 

Several  sample  cases  are  Included  on  the  following  pages  to  provide 
examples  of  the  output  for  the  several  options  available.  In  each 
case,  the  input  file  is  given,  followed  by  the  output.  The  sample 
cases  are  defined  in  Table  A-3. 


Table  A-3.  Sample  Cases 
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NAVAL  WEAPONS  CENTER,  CHINA  LAKE 


PROPERTIES  ALONG  RAT  TRAJECTORT 


trj  —  •~*-f\iryryojAj(Nj(\jrjrvjrgrj?Mr'jcgr\jcvjOqr\ir\jr\;rMr\jc\j<\jr\jr>jr'jrMfV  <\j  r*j  —  —  —  -- 

^  coooooooooooooocooooooooooooooooooooooooo 

X  444444444444444444444444444444444  +  4*44444 

_  >.i  f.i  t.i  ^.i  us  UJ  tjJ  u3  t%)  I«1  UJ  ur  i*J  U-*  l*-1  u>  (t)  1^,1  i«:  t.i  f.t  t.i  r.i  ;.i  i.»  r ■ ;  i».»  ui  ui  i.i  ;»i  u?  t.i 

— i  f7'^f\jOj.'\jcr-mr"iO  •-  mm<v'ir»*)Our'aDO  O'  la  -s  iAO'OEDir'Of^Tnwr,^o(,rirri7'i\H,'jM^C’ 

pof\nf>r»;3  wry  r^vO  ®  o  r-—  <\j  a  iT>0  Mio  n  Pn  r\j  la  rvm 

O  or-o  Oyo',/>0'C('jr~M^Ot\(rnr\iOP'ryf»f\jicOirip'SO'CnON^r^OO 


•-rvf'"UAc— ■’**-~CM<M<\if<na3-iAvA'©'Or— c-r-r-“f— OvOiALnaarAcurMoj*-  — 


x  00‘0'iAtA=Xt'"«C\i*-»-*“C\j<\irAacXJTlAiA'LAlAiAiAiAaoa  mCMAj*-  —  •-AjmaiAiAO'O'O 

^  — -oooooooooooooooooooooooooooooooooocoooo  — 

•-  I  I  I  I  t  I  I  I  I  1444444444444444444*44  1  I  I  II  III  II 

•  UJUJUJUJUiUJUJUJUJUJUJUJUJUJUJUJUJLtJUJUJUJUJUJUJUJUJUJUJUJUJUJUJUJUJUJUJUJUJUJUJU} 

«c  a  OOsO  f'OCT-rOr-f^tp'il'O  Ni£iO^SO(OONaiBCTNnO>ON'OOl^SS^CTP'''£0  O  a 

x  oOLnoiririr-coi^fMr^t:.  i\'Oc?'i/>oO’-cD”ncc-ooir'a'^r\jco^3-(\:a'CDP~irinO'jM» 

O  r-NCMr'(^ir»',nf^^v£)NO'0®^oo^pri3sco3  »  nc'{onco\coNom'ym'j\mlrip'Np- 


GO  ia*— f'-AjO'iA'-*— <\icD’-  c\j  ta  •--  ca  ca  <\j  ca  a  a  a  pa  r\j  o'  a-i  —  pa  \£>  Aj  *— aor\j  —  »— LAO'A.'r—'—  ir> 


o  oocooocoooooooooooooooooooooocoooooooooco 

X  4444444444444444444*4444  444  44444444444444 

<  UJ  111  UJ  UJ  UJ  UJ  UJ  UJ  UJ  UJ  UJ  UJ  UJ  UJ  UJ  UJ  UJ  UJ  UJ  UJ  UJ  UJ  UJ  UJ  UJ  UJ  UJ  UJ  UJ  UJ  UJ  UJ  UJ  UJ  UJ  UJ  UJ  UJ  UJ  UJ  UJ 

_j  p-ir>nff\r'in>D  os^oo  O'Ofvi  iMDfvi  f\jcD3rorv(Mttirr\jiccDo^'^NO'£Mrip-j'3  la  r— 

m  •5'-ff'f\.O0'O'S,-33a'P\jOc0ir'fVO®>O'CCC0OJ\ir'c0Cfy933'—  t—0'0'Oc\jO'»— a 

—  s-mr-u*"0r'*-ao0'0'0'0'0'0>cocoooee>r->t-t— r^f-eocococoo'O'O'O'O'O'CDr-vOiAt-p'Ar^  — 


t—  •—  pa  a  tAt/'iALAiAiAU''u‘'LAiAiAiALAkA,.AiAiA4AiAuAiAiAvAiALAiALAiALAiAiALA,.Aa  pa  p— 


X  AjrufyiCYCMiMrvj<MAjAj<MiM<M<\jr\j'\i<\ji\jfvjcMfM<\j<Mf'jr\jt\jAj(\:Ajc\jr\jrMr\j<\iAj<\ic\j;\irvj<\f\! 
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INTERSECTION  WITH  PLUME=  T 

ATTENUATION;  41.1665  DB 
PHASE;  6. ?8679E*03 
PHASE  DIFFERENCE:  -5. 898UHK-0 1 
DOPPLER  SHIFTED  FREOUF.NCT:  O.OOOE.OO  KHZ 


ANTENNA  GAIN  PATTERNS 
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VOLUME  INTEGRATED  SCATTERED  ENERGY  AT  DOPPLER-SHIFTED  FREQUENCY 
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FREQUENCY  INTEGRATED  SCATTERED  POWER  =  -163.660  DB 
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SCATTERED  AMPLITUDE, DB/HZ 


naval  weapons  center,  china  lake 

EXHAUST  PLUHE  RADAR  FREQUENCY  INTERFERENCE  CODE(PRFIC) 
WRITTEN  BY  AERONAUTICAL  RESEARCH  ASSOCIATES  OF  PRINCETON 
7-SEP-82 
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PROPER! I KS  ALONG  RAY  TRAJECTORf 
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INTERSECTION  WITH  PLUME-  T 

ATTENUATION:  60.6278  DB 

PHASE:  6.28666E.03 

PHASE  DIFFERENCE:  -6.93808E-01 
DOPPLER  SHIFTED  FREQUENCY:  P.000E*00  KHZ 


VOLUME  INTEGRATED  SCATTERED  ENERGY  AT  DOPPLER-SHIFTED  FREOUENCY 
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FREQUENCY  INTEGRATED  SCATTERED  P0WER=  -203-328  DB 
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SCATTERED  AMPLITUDE, DB/HZ 


NAVAL  WEAPONS  CENTER,  CHINA  LANE 
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PROPERTIES  ALONG  RAY  TRAJECTORY 
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ANTENNA  CAIN  PATTERNS 


VOLUME  INTEGRATED  SCATTERED  ENERGY  AT  DOPPLER-SHIFTED  FREQUENCY 
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FREQUENCY  INTEGRATED  SCATTERED  POWER:  -15N. 200  DB 
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SCATTERED  AMPLITUDE ( DB/ HZ 


TESTFF.DAT 

TEST  CASE  NUMBER  A  FOR  PRFIC,  MULTIPLE  CASES,  TRANSVERSE  ATTENUATION 
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PROPERTIES  ALONG  RAY  TRAJECTORY 
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PROPERTIES  ALONG  RAY  TRAJECTORY 
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Appendix  B 

SOFTWARE  DOCUMENTATION 


This  appendix  contains  the  software  documentation  for  the  Naval 
Weapons  Center  Plume  Radar  Frequency  Code  (PRFIC).  Brief  descriptions 
are  given  of  those  subprograms  that  contain  computations  of  physical 
quantities.  Those  subroutines  that  perform  input  and  output 
operations,  read  files,  or  provide  plots  are  not  described  in  detail. 
A  full  list  of  the  subroutines  was  given  in  Appendix  A  (Table  A-l). 
The  individual  descriptions  are  given  below,  followed  by  a  complete 
listing  of  the  entire  program. 

SUBROUTINE  PRFIAC( AANG, AGAIN, NA) 

Subroutine  PRFIAG  provides  the  antenna  gain  pattern  for  the 
transmitter  and  receiver.  AANG  is  the  angle  between  the  direction  of 
propagation  and  normal  to  the  antenna,  in  degrees.  AGAIN  is  the  gain 
at  this  angle,  in  dB,  and  NA  is  a  code  to  define  the  gain  distribution 
for  the  transmitter  (NA  ■  1)  and  receiver  (NA  *  2).  If  the  uniform 
distribution  option  is  chosen,  then 


AGAIN  -  0;  0  <  AANG  <  90. 


(B-l) 


If  the  gain  distribution  is  nonuniform,  a  linear  interpolation  between 
angles  is  used.  If  G^  is  the  gain  at  the  angle  9^,  then 


AGAIN  -  Gt  -  (9i  -  AANG)  x  (Gi  -  Gl+j)/(9i  -  9i+1)  (B-2) 

where  the  input  order  is  such  that  9^+j  >  9^.  If  AANG  is  beyond  the 
range  of  the  last  angular  entry,  the  subroutine  returns  the  gain  at  the 
last  entry. 


SUBROUTINE  PRFIDT( ATTEN, PHASE , DOPPS) 

PRFIDT  is  the  primary  subroutine  that  governs  the  direct  transfer 
of  radiation  from  the  transmitter  to  the  receiver  along  the 
line-of-sight  of  a  refracted  ray  trajectory.  It  returns  the  total 
attenuation  (ATTEN)  and  phase  (PHASE)  for  the  path  and  the  Doppler 
shift  of  this  transmitted  ray  (DOPPS). 
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The  primary  function  of  PRFIDT  is  to  set  up  the  paths  exterior  and 
interior  to  the  plume  and  then  to  call  PRFIIT  to  do  the  integration  for 
each  segment.  However,  the  first  task  is  to  determine  if  direct 
transfer  can  occur.  No  direct  transfer  is  allowed  for  the  three  cases: 


1.  co-located  transmitter  and  receiver  (monostatic  radar) 

2.  angle  between  the  normal  to  the  transmitter  and  the  direction  to 
the  receiver  >  90  degrees 

3.  angle  between  the  normal  to  the  receiver  and  the  direction  from  the 
transmitter  >  90  degrees 


If  any  of  the  conditions  occur,  the  subprogram  returns  and  prints  a 
message. 

If  the  straight  line-of-sight  option  is  used,  the  direct  transfer 
occurs  in  three  steps:  (1)  transmitter  to  plume;  (2)  through  the 
plume;  (3)  plume  to  receiver.  PRFIPI  is  used  to  define  the 
intersection  with  the  plume  (one  intersection  or  two,  if  the 
line-of-sight  enters  the  plume).  If  there  is  no  intersection,  the 
program  returns  the  attenuation  (»  0),  phase  and  Doppler  shift,  and 
prints  a  flag,  INTERSECTION  WITH  PLUME  =*  F. 

If  the  ray  trace  option  is  used,  an  iterative  solution  is  required. 
The  initial  guess  for  the  curved  trajectory  is  taken  to  be  the  straight 
line-of-sight.  At  the  endpoint  of  the  integration,  an  error  vector 
between  the  computed  endpoint  and  actual  endpoint  (receiver  location) 
is  calculated.  If  this  error  exceeds  a  specified  criterion  (EPSA(l)), 
then  a  new  guess  for  the  initial  direction  at  the  transmitter  is  made. 
The  current  new  guess  for  the  direction  cosines  at  the  transmitter  is 


(cosa)n+1  -  (cosa)n  +  RC  *  ERV(I)/SL  (B-3) 


where  n  is  the  old  guess,  RC  is  an  overrelaxation  factor 
(0.5  <  RC  <  1.0),  ERV(I)  is  the  error  vector  from  the  computed  endpoint 
to  the  desired  endpoint,  and  SL  is  the  length  of  the  ray  path.  Similar 
expressions  apply  for  the  other  two  angles  g(I  »  2)  and  -y(I  »  3).  This 
is  a  simple  iterative  solution  to  the  two-point  boundary  value  problem 
that  adjusts  the  initial  ray  direction  to  achieve  convergence.  The 
iteration  is  performed  until  convergence  occurs  or  until  the  number  of 
iterations  exceeds  NIT,  an  input  variable.  If  no  convergence  occurs,  a 
warning  message  is  printed. 
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SUBROUTINE  PRFIEP 

The  subroutine  PRFIEP  utilizes  the  input  flowfield  variables  read 
directly  from  the  SPF  input  file  to  create  an  array  of  quantities  used 
in  the  RF  interference  calculations.  It  is  called  once  at  the 
beginning  of  each  run  when  a  new  flowfield  file  is  defined.  The  RF 
interference  properties  are  then  calculated  at  each  point  in  the 
flowfield,  once  and  for  all.  All  subsequent  interpolations  required  in 
the  attenuation  and  scattering  calculations  throughout  the  plume  volume 
are  performed  from  these  quantities  directly,  rather  than  for  the 
flowfield  variables  initially  on  the  input  file.  The  objective  of  this 
calculation  sequence  is  to  minimize  the  number  of  operations.  For 
example,  calculating  the  RF  properties  first,  and  then  interpolating, 
requires  fewer  calculations  than  interpolating  for  the  flowfield 
properties  and  then  computing  the  RF  properties  at  each  point  along  a 
ray  trajectory. 

The  input  flowfield  file  is  stored  in  the  array  FF(I,J,K).  PRFIEP 
calculates  the  RF  interference  properties  and  stores  them  back  in 
FF(I,J,K).  The  individual  quantities  are  defined  below. 


Input  Flowfield  properties 

FF(I.J.K) 

Units 

1 

r , 

radius  at  axial  location  z(I),  radial  location  J 

m 

2 

T, 

temperature 

K 

3 

P, 

pressure 

n/ 

4  ne ,  mean  electron  number  density 

5  ng^/n| ,  normalized  mean  square  electron  number 
density  fluctuation 

6  ven,  electron-neutral  collision  frequency 

7  A,  turbulence  macroscale 

8  u,  mean  axial  velocity  m/s 

9  q2/u2,  normalized  mean  square  velocity  fluctuation  - 


|C 

Aa 
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This  is  a  multiplier  that  is  common  to  the  scattering  cross  section 
used  for  both  the  attenuation  and  scattering  calculation.  For 
scattering,  o  multiplies  the  angular  dependent  scattering  function.* 


_ sin2x 

[l  +  4(kA)2sin*  !• 


z  e/Ul/S 
2 


( B-5) 


For  extinction,  a  multiplies  the  integral  of  this  function  over  all 
directions.  This  integral  is  cal  Led  Q  and  depends  on  kA.  It  Is 
evaluated  in  function  subprogram  SCAT(kA). 

K  *  7,8  u  and  q2/u2  carried  over  directly  from  flowfield. 

K  *  9,10  Radial  and  axial  components  of  the  index  of  refraction 

gradient . 

The  index  of  refraction  gradients  for  use  with  the  ray  tracing 
option  are  calculated  according  to  the  following  finite  difference 
approximations.  Subscripts  i,j  denote  z,r  locations,  respectively. 


(—) 

,  1  ' 

Ur/i.j 

2 

(B-6) 


'3*/ i,j 


1  i+l’ 

2 


zi+l  "  zi 


+  n±i)  ~  "i-l.J 

zi  '  zi-l 


( B-7) 


SUBROUTINE  PRF1 1T( X0 , SO , XF , ATTEN , PHASE, ERV , 
ERVM  ,XFV,  SFV,  IN',  RAY) 


Subroutine  PRFi.IT  computes  the  ray  path  between  two  points  either 
outside  or  inside  the  plume.  If  the  ray  trace  option  is  chosen,  the 
ray  trajectory  within  the  plume  is  computed  with  the  geometric  optics 
approximation.  If  a  line-of-sight  option  is  chosen,  the  ray  trajectory 
is  a  straight  line.  In  both  cases,  the  integrated  attenuation  and 
phase  along  the  path  are  computed. 


*For  one  choice  of  the  turbulence  spatial  correlation. 
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XO(I),  1-1,3  :  Coordinates  of  initial  point,  m 
S0(I),  1-1,3  :  Direction  cosines  of  initial  direction 


XF( I) ,  1-1,3 


ATTEN 


PHASE 


Coordinates  of  final  point  of  the 
line-of-sight 

Integrated  attenuation  along  the  path 
Integrated  phase  along  the  path 


ERV(I),  1=1,3:  Error  vector,  ERV(I)  -  XF(I)  -  XFV(I) 


straignt 


Magnitude  of  error,  ERVM  -  I  T*1  ERV(  I) 


XFV(I),  1*1,3:  Final  point  of  integration 


SFV(I),  1-1,3:  Direction  of  path  at  the  endpoint 


IN  -  T,  trajectory  within  the  plume 
»  F,  trajectory  outside  the  plume 

RAY  *  T,  ray  trace  calculation 

-  F,  line-of-sight  calculation 


Exterior  to  the  plume,  the  attenuation  and  phase  are  given  simply 
by  ATTEN  =  0,  and 


PHASE  -  2ttv0s/c 


(  B-8) 


where  s  is  the  path  length  between  the  initial  and  final  point.  The 
endpoint  error,  ERVM,  is  zero,  the  final  points  are  equal, 
XF( I)  =  XFV(I),  and  the  final  direction  and  initial  directions  are  the 
same  S0(I)  -  SFV(I).  This  solution  applies  regardless  of  the  ray  trace 
option. 

Interior  to  the  plume,  the  line-of-sight  is  divided  into  increments 
according  to 


AS  -  s/NCUTS 


(  B-9) 


If  the  line-of-sight  option  is  taken,  the  integration  for  the 
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attenuation  and  phase  is  summed  along  the  straight  1 ine-of-sight  by  a 
trapezoidal  rule 


NCUTS 

ATTEN  -  £  T(“i+1  +  ai)  (B-10) 

i-1  2 


NCUTS 

PHASE-  £  T(Si+l  +  8i}  C  B— 11) 

i-1  2 


In  both  cases,  the  attenuation  and  phase  constants  are  obtained  at  each 
point  i  by  interpolation  between  the  values  stored  in  the  flowfield 
file. 


If  the  ray  trace  option  is  chosen,  then  the  integration  within  the 
exhaust  plume  proceeds  along  s,  the  arc  length,  but  this  no  longer 
coincides  with  the  straight  line-of-sight .  The  equations  defining  the 
ray  trajectory 

r(s)  —  | x ( s )  ,  y(s) ,  z(s)}  ( B- 1 2 ) 

x(s  +  As)  »  x(s)  +  As  +  -M~-  -  ~~(f  •  Vn))  (As)2 

vds  /  2n  '3x  ds  ' 

8  s  ( B-l  3) 


y(s  +  as) 


y(s) 


As  + 


wjm. 

2n  \  3y 


Vn))  (As)'1 

(B-14) 


z(s  +  As) 


z(s) 


As  + 


_l/3n 
"in'  3z 


pit 

ds 


Vn)j  (As)2 

s  (B-l 5) 


where 


dx  cTy  dz 
d  s  ’  d  s  ’  cTs 


and 


(B-16) 
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Vn 


cos8 


3n 


,  a  3n 

sinQ  — , 
3  r 


3n 

3z 


(B-17) 


Solution  of  Eqs.  ( B— 1 3 ) — ( B— 15)  is  achieved  in  a  stepwise  numerical 
integration  throughout  the  plume.  The  initial  point 


r(0)  -  x( 0) ,  y(0) ,  z(0)  (B-18) 


and  direction 


t(0)  -  ^  (0).  ^  (0).  j-  (0)  ( B— 19) 

ds  ds  ds 


are  known.  However,  for  a  given  initial  point  and  direction,  the 
endpoint  and  direction  are  unknown  until  the  integration  along  the 
trajectory  has  been  performed.  For  a  .  nown  initial  point  and  known 
endpoint,  the  solution  for  the  ray  trajectory  is  a  two-point  boundary 
value  problem.  Ai  iterative  solution  must  be  performed  until  the 
initial  direction  is  found  which  allows  the  trajectory  to  pass  between 
these  points. 

This  subroutine  does  not  perform  the  iteration  but  simply 
integrates  the  ray  equations  between  two  points.  The  actual  choice  of 
a  new  direction  is  made  in  the  individual  subroutines  PRFIDT  and  PRFISC 
(for  the  directly  transmitted  ray  and  scattered  ray,  respectively). 
Integration  along  the  trajectory  proceeds  until  the  distance  between 
the  known  endpoint  XF(I)  and  point  on  the  ray  path  is  smaller  than 
As/2.  A  check  is  included  to  insure  that  if  the  trajectory  does  not 
pass  within  As/2  of  the  known  endpoint,  the  integration  stops  when  the 
error  is  least.  When  this  criterion  is  met,  the  subroutine  returns  the 
actual  endpoint  of  the  integration  XFV(I),  the  actual  direction  of  the 
ray  at  the  endpoint,  the  error  vector  ERV(I),  and  magnitude  of  the 
displacement  error  ERVM.  These  quantities  are  used  to  choose  a  new 
initial  direction.  If  the  line-of-sight  option  is  chosen,  the  errors 
are  zero  and  the  final  direction  is  the  same  as  the  initial  direction. 


SUBROUTINE  PRFILI( Z , R, VZR) 

The  subroutine  PRFILI  is  a  linear  interpolation  program  that 
computes  the  RF  interference  properties  at  the  axial  location,  z,  and 
radial  location,  r  which  do  not  coincide  with  the  grid  points  of  the 
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input  flowfield.  Let  i,J  be  indices  identifying  the  axial  and  radial 
locations,  respectively.  Then,  if  f  is  any  dependent  variable  such 
that  f^j  is  its  value  at  the  grid  points  i , j  and  f2>r,  the  value  at 
(z,r): 


‘-•r '  '  '“‘•J-) 


( B-20) 


r  -  r 


£i,r  ‘  fi, 


+  ‘  (f  .  -  f  \ 

,j_I  ri,j  “  rt,j-l  ■  i,J  l,j~1' 


( B-2 1 ) 


then 


z  -  z 


f„  _  -  f,  .  _  + 

2,r  i— l,r  2  —  2 


i- 1 


‘i  “i-1 


(<l,r  -  <>-22> 


The  values  of  the  nine  quantities  n,  a,  6,  kA,  a,  u,  q2/u2,  3n/3r, 
3n/3z  are  returned  in  the  array  VZR  for  each  point  z,r,  even  though  not 
all  are  used  at  each  call.  For  example,  3n/3r  and  3n/3z  are  not  used 
when  the  1 ine-of-sight  option  is  chosen. 

If  the  point  is  outside  the  plume  boundary,  ambient  values  are 
returned.  The  non-zero  quantities  are  n-1,  |}=2nv/c,  and  u=u<„. 


SUBROUTINE  PRFIPC(v,  vp,  ven,  a,  S,  n) 

fhe  subroutine  PRFIPC  calculates  the  plasma  index  of  refraction, 
attenuation,  and  phase  constants,  given  the  radar  frequency,  v,  plasma 
frequency.  y_,  and  electron-neutral  collision  frequency,  ven. 


vp  -  8.97  n^2  ( B-23) 

wp  -  2nVp  (B-24) 

(i)  *  2nv  (  B-2  5) 
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iKj  »  1  - 

(iDp/w)2 

_  ^  <* 

<« >>2  ^ 

*  U> 

( B— 26 ) 

\  / 

“  T  1 

I2 

■ +  (^)2 

|K|  - 

(kr  +  Ki>1/2 

(B-27) 

-I'M 

II 

c 

|k|  +  kr) 

1/2 

( B-28) 

<D  1 

a  * - 

c  2 

(|x|  -  Kr)! 

J  1/2 

( B-29) 

0 

=*  nu)/c 

( B-30) 

SUBROUTINE  PRFIPI( XO , XF , XI , XE , MTYPE , NTYPE) 

It  is  convenient  to  separate  the  portions  of  the  propagation  path 
that  are  outside  and  inside  the  plume.  Subroutine  PRFIPI  solves  for 
the  points  at  which  a  straight  line-of-sight  intersects  the  plume 
boundary.  This  subroutine  is  used  to  answer  the  following  questions: 


1.  Does  the  line-of-sight  intersect  the  plume? 

2.  If  there  is  an  intersection,  how  much  of  the  line-of-sight  is 
outside  the  plume?  The  attenuation  and  phase  of  the  outside 
portion  can  be  computed  in  one  step  because  it  is  a 
homogeneous  path  (with  a  -  0  and  0  *  2irv/c),  while  only  the 
nonhomogeneous  portion  interior  to  the  plume  need  be  broken 
into  smaller  steps  for  integration. 

a.  This  calculation  completely  defines  the  propagation  path 
when  refraction  is  not  considered. 

b.  This  calculation  gives  a  first  guess  for  the  ray 
trajectory  when  refraction  is  considered. 
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The  solution  is  determined  by  the  following  sequence  of  operations. 
Let  =*  (xj.yj.zj)  be  the  initial  point  and  ^2  *  (x2»y2>z2^  be  the 
final  point  along  the  path.  Then  any  point  along  the  path  between 
and  is  t  »  (x ,y  ,z) , 


r  »  Xj  +  ^  s 


(B-31) 


where  5  is  a  scalar  distance  (£  »  0)  and  s  is  the  unit  vector  in  the 
direction  of  t 2  from 


(r. 


-  V' 


-  ?, 


(B-32) 


The  radius  of  the  plume  at  any  axial  location  is 


r(z)  -  rt 


+ 


ri+l  -  ri 
zi+l  "  zi 


(z  -  zi ) 


zi  <  2  <  zl+1 


(B-33) 


where  r^.z^  are  the  plume  radius  and  axial  location  at  the  ich  station 
given  in  the  flowfiel  1  file  li.e.,  FF(  I , NPTS( J) , 1 ) ,  Z(I),  1  <  I  <  NZ ) . 
At  the  point  of  intersection 


(? 


-  2 
i) 


+  (r 


(xj  +  £s 


i)  +  ( y  j  +  $s 


*  2 
j) 


r2(z) 


( B-3A) 


and 


?•  k-Zj+£s  •  k  »  z  ( B—  3  5) 

where  I ,  j  ,  ft  are  unit  vectors  in  the  x,  y,  and  z  directions.  These 
two  equations  can  be  combined  to  give  a  single,  quadratic  equation  for 


af,2  +  2b£  +  c  -  0 


( B-36 ) 


with 
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a  «  B2  - 

(s  • 

i)2  -  (s  .  j)2; 

ri+l  "  ri 

A  ,  r  +^±! - ±  (*  -  Zl) 

zi+l  "  zi 

(B-37) 

b  -  AB  - 

X.  (s 

•  i)  -  (s  •  j); 

b  -  ri+1  '  ri  (;  .  w 

( B-38) 

1 

zi+i  -  zi 

c  * 

A2  -  (Xj  +  yj) 

( B-39) 

and 

£  -  (-b  +  /b2  -  a  c  ) / a  (B-40) 

The  two  roots  to  this  equation  are  evaluated  for  each  conical  segment 
of  the  plume,  T  »  1 . NZ.  The  solution  of  interest  must  satisfy 


Zi  «  2  <  zi+l  ‘  (B"41) 

That  is,  the  line-of-sight  must  pass  through  the  conical  segment  of  the 
plume  at  which  the  radius  is  defined  by  Eq .  (B-33).* 

The  exceptional  cases  occur  when 


b2  -  ;:c  <  0 


(tf-42) 


for  which  there  is  no  intersection,  and  when  a  =*  0,  for  which  there  is 
only  one  solution  (the  line-of-sight  can  intersect  the  plume,  but  the 
angle  of  the  plume  expansion  is  such  that  it  will  never  exit).  In  the 
event  that  the  line-of-sight  does  not  exit  the  plume,  we  define  the 
line-of-sight  interior  to  the  plume  to  extend  to  the  end  of  the  input 
flowfield.  (This  situation  can  occur  if  the  receiving  antenna  is 
downstream  of  the  end  of  the  flowfield  and  sufficiently  close  to  the 
plume  axis.) 


*It  is  possible  for  the  line-of-sight  to  intersect  the  upstream  and 
downstream  extensions  of  the  conical  volume  defined  by  Eq .  (B-33). 
That  is,  5  is  a  root  of  Eq.  (B-40),  but  z  <  z^  or  z  >  z^+j.  However, 
that  is  not  a  valid  solution  for  the  plume  intersection. 
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In  general,  there  are  two  real  roots  that  result  from  the  solution 
after  all  plume  segments  are  examined.  From  the  definition  of  6  *  0, 
the  root  defining  the  incident  intercept,  is  the  smaller  of  the  two, 
Ci  <  £e,  where  £e  defines  the  exit  point.  The  positions  of  the 
intercepts  are 


X1  +  ?i,e 


( B-43) 


The  same  formulation  applies  if  the  initial  point  is  outside  the 
plume  and  the  endpoint  is  inside,  or  vice  versa.  This  situation  occurs 
in  the  volume  scattering  calculation.  If  is  outside  and  *2  inside, 
then  the  single  root  of  interest  is  ^ .  If  is  inside  and  *2 
outside,  then  one  of  the  two  roots  is  negative,  and  we  take  the 
positive  root  to  be  £e. 

All  these  cases  are  computed  in  the  subroutine  PRFIPI,  and  are 
defined  in  the  arguments,  as  follows: 


X0  ■  X0(I),  I  *  1,3  x,y,z  coordinates  of  initial  point 

XF  »  XF(  I ) ,  I  —  1,3  MTYPE  •'  0;  XF  *  direction  cosines  of  direction  to 

endpoint 


MTYPE  >  0;  XF  *  coordinates  of  endpoint 
XI  *  XI  (  I )  ,  I  -  1,3  coordinates  of  enc-y  point 

XE  XE(1),  I  *  1,3  coordinates  of  exit  point 

MTYPE  *  flag  for  type  of  input  (coordinates  or 

direction  to  final  point) 


NTYPE  *  flag  for  type  of  path 

NTYPE  -  1  Initial  and  endpoints  outside  the 

plume 

(solutions  are  XI  and  XE,  or  no  intersection) 
NTYPE  »  2  Initial  point  outside,  final  point 

inside,  solution  is  XI. 

NTYPE  «  3  Initial  point  inside,  final  point 

outside,  solution  is  XE. 
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SUBROUTINE  PRFIRZ( ZP, RP) 

The  integration  in  the  radial  direction  of  the  volume  scattering 
integral  (see  description  of  subroutine  PRFISC)  uses  a  radial 
coordinate  normalized  by  the  local  plume  radius  r(z).  Subroutine 
PRFIRZ  interpolates  between  neighboring  flowfield  stations  z^,  z^+i  to 
define  the  plume  radius  at  the  axial  quadrature  point  ; 
zl  <  zi  4  zl+l‘  A  simple  linear  interpolation  is  used: 


r(zi)  »  r(zj)  + — - - — fr(z4+1)  -  r(z4)l  (B-44) 

zl+ 1  "  zl  L  J 

where  r(zj)  =»  FF[ L, NPTS(  L.) ,  1 1 ,  the  outermost  radial  location  defined  in 
the  flowfield  file  at  the  axial  location  z(L). 


SUBROUTINE  PRFISC 

Subroutine  PRFISC  computes  the  scattered  signal  due  to  fluctuations 
in  the  electron  number  density.  It  performs  a  numerical  evaluation  of 
the  scattering  from  the  entire  plume  volume.  This  requires  a 
three-dimensional  volume  integration;  the  plume  volume  is 
axisymmetrlc ,  but  the  arbitrary  locations  of  transmitter  and  receiver 
make  the  integrand  dependent  upon  all  three  coordinates  (r,9,z). 

The  approach  used  in  this  report  is  to  evaluate  the  volume  integral 
by  a  succession  of  Gaussian  quadratures.  For  a  fixed  number  of  points, 
these  schemes  yield  the  most  accurate  quadratures.  The  volume  integral 
to  be  evaluated  is 


C(r ,9  ,z,v)  dV  (B-45) 


where  P(v)  is  the  frequency  spectrum  of  the  mean  square  scattered  power 
and  C  is  the  scattering  cross  section  at  the  point  r.9,z.  The 
integrand  also  depends  on  the  scattering  angle  between  the  directions 
of  the  incident  and  scattered  waves  and  also  the  distances  from  the 
scattering  location  to  both  the  transmitter  and  receiver.  In  the  event 
of  refraction,  a  separate  curved  ray  trajectory  must  be  computed 
between  the  point  r,9,z  and  both  the  transmitter  and  receiver.  In 
addition,  the  attenuation  and  phase  for  each  path,  for  both  the  curved 
ray  trajectory  or  straight  line-of-sight ,  is  included.  The  detailed 
specification  of  C(r,9,z,v)  is 
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C(r 


e  z  ,  A2  CTCA°  sinX  PDF(v  -  v0) 

8"  rV  .  ^  1  •  2  6lU/t 

rTrR  1  +  4(kA)  sin  — 

2 


(B-4b) 


where  Gj,  CA,  and  rT,  rR  are  the  antenna  gains  and  distances  from  the 
sea  te* *ng  point  to  the  transmitter  and  receiver,  respectively* 

In  FRFISC,  the  volume  integral  is  evaluated  according  to 

/Zf  /*R(z)  SI* 

I  I  C(r,9  ,z)rd9  dz  (B-47) 

i.i  */0  Jn 


E  E  E  C(rj»Vzi)rj  wk  wj  ui  (B-48) 

i  J  k, 


(zf  2i.)  »  J2  r2<zi>  E  nj  c(ni»Sk’H)wk  wi  wi 
i-I  j  J  k-1  J  k  j  1 


( B-4  9) 


^®re,  *i  *  WJ*.wkare  the  weights  and  Ck,  „ , ,  Ai  are  the  roots  of 
normalized  variables  in  the  9,  r,  and  z  directions,  respectivelv.  The 
normalization,  ranges,  and  number  of  points  are  given  in  the  following 
ta  e.  In  all  three  cases,  the  independent  variables  are  scaled  to  the 
interval  (-1,1)  so  the  quadrature  points  are  roots  of  Legendre 
polynomials.  & 


Table  B-i.  Numerical  Quadratures  Used  in  Volume 
Scattering  Calculation 


- - - - 

1 

;  Coordinate 

Transformed 

- - - 

Range 

9 

9  -  n(l-;) 

-1  t  4  <  1 

r 

r  «  n  R(z) 

0  t  q  <  1 

(1/2  of 

z 

z  -  iUf-zt) 

0  4  £  4  1 

- - - 

(1/2  of 

No.  of  Points 


10 


10 


24 
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The  axial  v  riable  is  scaled  by  an  initial  and  final  plume  length  , 
Zf,  respectively,  that  can  be  specified  in  the  input.  This  allows  the 
computation  of  scattering  from  a  finite  length  of  plume  different  from 
the  length  of  flowfield  in  the  input  file. 

The  integral,  Kq.  (B-47),  is  also  the  temporal  power  spectral 
density  of  the  scattered  signal,  i.e.,  P  =  P(v),  that  accounts  for  the 
probability  of  both  the  spatial  and  temporal  fluctuations  in  electron 
number  density.  The  integrand  contains  a  frequency  dependent 
multiplier,  the  probability  density  function  PDF(v).  p(v)  is  evaluated 
at  increments  in  frequency  difference,  Av  =  v  -  v0,  from  the 
transmitter  frequency,  v0 . 

If  the  ray  trace  option  is  chosen,  the  incident  and  scattered  ray 
paths  are  curved.  Iterative  solutions  are  applied  to  both  the  segments 
in  PRFISC.  The  procedure  is  like  that  outlined  for  the  directly 
transmitted  ray  in  PRFIDT. 


FUNCTION  PDF  (OF  ,VMK,QQ ,  ANG ,  ISFRUQ) 

Function  subroutine  PDF  returns  the  probability  density  function 
for  the  Doppler  shifted  scattered  wave.  The  Doppler  shift  includes  the 
effect  of  the  resultant  mean  velocity  of  the  transmitter,  nissile,  and 
exhaust  plume  relative  to  the  receiver.  There  is  a  frequency  spread 
about  this  mean  shift  due  to  the  turbulent  fluctuations.  At  present, 
it  is  assumed  that  the  velocity  fluctuations  are  isotropic  with  respect 
to  direction  at  each  point  and  have  a  Gaussian  distribution  with  a  mean 
square  value  of  q ' 2  about  the  local  time  averaged  velocity.  With  these 
assumptions,  the  probability  density  function  of  the  Doppler  shifted 
scattered  wave  at  a  point  is 


PDF( v  -  vQ) 


1 


2  ^[Zsinfo/Z)]2  2v0sin(9/2) 


( B-50) 


DF  «  frequency  difference  (v  -  vQ),  Hz 

VMK  *  component  of  resultant  velocity  of  transmitter,  missile,  and 
mean  plume  flow  relative  to  the  receiver,  m/s. 

QQ  ■  mean  square  turbulence  intensity,  q '  2  ,  m2/s2. 


ANG  »  magnitude  of  the  resultant  scattering  wave  vector,  2sin(0/2). 
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BFREn  »  transmitter  frequency,  v0. 

PDF  =  probability  density  of  mean  square  scattered  power  at  the 
frequency,  v ,  0  <  PDF  <  1. 

SUBPROGRAM  SCAT(kA) 

The  function  subroutine  SCAT(kA)  calculates  the  directionally 
averaged  scattering  component  of  the  electron  number  density 
fluctuation.  It  depends  on  the  single  augment  kA  and  is  given  by 


SCAT  -  ~\^  (2p2  +  2P  +  l)|l  -  (1  +  2p)_5/b  J 


2(p  +!)(!+  2p)1/t3  -  l|  +  yj(l  +  2p) 7/b  -  ijj  ( B-51 ; 


where  p  =>  2(kA)^. 


FUNCTION  SCAI.P( A , B) 


Function  subroutine  SCALP  returns  the  scalar  product  of  two  vectors 
defined  in  rectangular  Cartesian  coordinates. 


3 

SCALP  -  £  A(I)  *  B(  I )  ( B-52) 

I»1 


A(I)  and  B(I),  where  I  »  1,2,3,  are  the  x,  y,  and  z  components  of  the 
two  vectors,  respectively. 
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Code  Listing 


[ 
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S  BROUTINE  PRFIAC(AANG, AGAIN, NA) 

C 

C  COMPUTES  ANTENNA  GAIN  AT  ANGLE  AANG 

C 

C  AANG-ANGLE ,  DEG 

C  A'AIN-GAIN,  DB 

C  SA»1,  TRANSMITTER 

C  *2,  RECEIVER 

C 

I NCLUDE ' PRF1CM. FOR/  NOLI  ST’ 

Ir(NANC(NA).NE.O)  GO  TO  5 
AGATN*0. 

RETURN 

5  C  NTINUE 
N>NANG(SA)-1 
D.  10  1*1, NC 

IF  AANG. GE. ANGLE (NA, I ) .AND. AANG. LT. ANGLE!  NA, I+l ) )  GO  TO  20 
10  C  NT I SUE 

AGA IN*CA I N( NA , NANC( NA ) ) 

RETURN 

20  C  NTINUE 

AGA IN*GA  IN(NA,I)-(  ANGLE!  NA , I ) -A ANG ) * ( GA I N( NA , I )  -GA  IN( NA , 1+ 1 ) ) / 
1  ANCLE ( NA , I ) -ANCLE ( NA , l+l)) 

RETURN 


*********************************************************************** 

C  PRrlCM 

r' 

V* 

C  «:  MMON  FOR  PROGRAM  PKFIG 

C  FLUME  RADAR  FREQUENCY  INTERFERENCE  CODE 

C  WRITTEN  FOR  NWC,  CHINA  LAKE 

C 

LOGICAL  RAYTR , FSCAT , INSECT, PATH, EXTER 
REAL  NCUTS 
BYTE  FILEFF( 20) 


C 


PARAMETER  ( PI*3 . 14 i 592653 ,  C-2.998E+08) 

D'  NOT  ALTER  THE  POSITIONS  OF  VARIABLES  IN  THE  NEXT  8  LINES 

C ’ MMON  PFREQ.ZLL.ZUL 
CIMMON  XT( 3) ,DCT( 3) ,DCT0( 3) 

C " MMON  VT( 3) ,XR( 3) ,DCR( 3) 

C  j'TMON  DCRO(  3) ,  VR(  3)  ,XM(  3) 
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COMMON  DCM( 3) ,DCMO( 3) ,VM( 3) 
COMMON  NCUTS,RC,NIT 
COMMON  EPSA(2) ,DUMX 
COMMON  RAYTR , FSCAT 

C  OK  ... 


COMMON  RTIIOM ,  INSECT ,  PATH ,  EXTER 

COMMON  TITLE( IS) .FILEFF 

COMMON  FF(60, 50,  10)  , 7.(60)  ,NPTS(60) 

COMMON  VZRAMB( 9) , ILAST,NZ 

COMMON  NANG(  2)  ,ANC.I.E(2 , 50)  ,GAIN(  2 , 50) 

C  END 

-k  ■kit -kit  iekicyfk'kick’k’k'k^ifk-k  "k-kifkiclrkic-k-k  ★★★****  *****★★*  *  **★*★★*★★***★★*★★★★★ 

SUBROUT I NE  PRF IDT( ATTEN , PHASE , DOPPS ) 

C 

C  CALCULATE  RAY  PATH  FROM  TRANSMITTER  TO  RECEIVER 

C 

INCLUDE  ’ PRF I CM. FOR/ NOLI ST* 

LOGICAL  IN, RAY 

DIMENSION  SU( 3) ,DUM( 3) ,ERV(3) ,XFV( 3) ,SFV( 3) ,XD(3) , X I ( 3 ) , 
i  XE(  3) ,XF(  3) , VTR( 3) 

C 

DETERMINE  IF  DIRECT  TRANSFER  CAN  OCCUR 


PATH=.TRUE. 

DOPPS=0. 

SS=0. 

DO  10  1=1,3 

SS=SS+(XR( I)-XT( I))**2 
10  CONTINUE 
SL=SQRT( SS) 

ATTEN=0. 

PHASE=2.*PI*SL*BFREQ/C 

IF( SS.GT.O. )  GO  TO  30 
WRITE(6 , 20) 

20  FORMATOH  COLOCATED  ANTENNAS') 
INSECT=. FALSE. 

RETURN 

30  CONTINUE 
DO  40  1=1,3 

SU(I)=(XR(I)-XT(I))/SL 
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VTR( 1)-VT(L)-VR(  I) 

40  CONTINUE 

COSPT-SCALP!SU,DCT) 

LF!COSPT.GT.O. )  GO  TO  60 
WRITE(b,50) 

50  FORMAT! 1H  , ’ NO  DIRECT  PATH  TO  RECEIVER,  S*NT<0. * ) 
INSECT-. FALSE. 

RETURN 
60  CONTINUE 

CuSPR— SCALP!  SU ,  OCR) 

1F!COSPR.GT.O. )  GO  TO  80 
WR1TE!6,70) 

70  FORMAT! 1H  , ’ NO  DIRECT  PATH  TO  RECEIVER,  S*NR<0.') 
INSECT-. FALSE. 

RETURN 
80  CONTINUE 

C  DOPPLER  SHIFT  OF  DIRECT  TRANSMISSION 

VTMS-SCALP! VTR , SU ) 

DOPPS-BFREQ*VTMS/C 


IF! PATH)  WRITE(6 , 300)  !XT(I) ,1-1 ,3) 

C 

C  INTERSECTION  WITH  PLUME 

C 

C  LOS  THROUGH  PLUME 

C 

90  CONTINUE 

IF(RAYTR)  GO  TO  100 

CALL  PRFIPl! XT ,XK ,XI ,XE, 1 , 1 ) 

IF( INSECT)  THEN 
PHASE-0. 

EXTER-.TRUE. 

IN-. FALSE. 

RAY-. FALSE. 

CALL  PRF 1 IT! XT , SU , XI , ATTENT , PHASET , ERV , ERVM , XF V , SFV , I N , RAY ) 

ATTEN-ATTEN+ATTENT 

PKASE-PHASE+PHASET 

C 

C  LOS  CALCULATION 
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C 


IN-. TRUE. 

CALL  PRFI IT( XI , SU , XE , ATTENT , PHASET , ERV , ERVM , XFV , SFV , IN , RAY ) 

ATTEN-ATTEN+ATTENT 

PHASE=PHASE+PHASET 

EXTER=. FALSE. 

I N=. FALSE. 

CALL  PRFIIT(XFV,SU,XR, ATTENT, PHASET, ERV, ERVM, XFV, SFV, IN, RAY) 

ATTEN=ATTEN+ATTENT 

PHASE=PHASE+PHASET 

DB=ATTEN*8. 6858 

I F ( PATH )  WRITE(6 , 30U)  (XR( I) , 1-1 ,3) ,  DB, PHASE 

END  IF 
■RETURN 

RAY  TRACE  CALCULATION 


100  CONTINUE 
IT-1 

PATH-. FALSE. 
105  CONTINUE 


ATTEN-O. 

PHASE-0. 

OUTSIDE  PLUME 

EXTER-.TRUE. 

IN-. FALSE. 

RAY-. FALSE. 

CALL  PRFI PI (XT, SU,XI,XE,-1 , 1) 

CALL  PRF 1 IT( XT , SU , XI , ATTENT , PHASET , ERV , ERVM , XFV , SFV , IN , RAY ) 

ATTEN- ATTE N+ATTE NT 

PHASE-PHASE+PHASET 

INSIDE  PLUME 


IN-. TRUE. 

RAY-. TRUE. 

CALL  PRF I IT( X I , SU , XE , ATTENT , PHASET , ERV , ERVM , XFV , SFV , I N , RAY ) 
ATTEN-ATTEN+ATTENT 
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C 

C  OUTSIDE  PLUME 
C 

SA=0. 

DO  110  1=1,3 
SA=SA+(XFV( l)-XR( l))**2 
DUM( I)=SU( I) 

110  CONTINUE 

SA=SQRT(SA) 

DO  120  1=1,3 
XF(I)=XFV( I)+SA*SFV( I) 

XD( I )=XFV( l) 

SL'(  I  )=SFV(  I ) 

120  CONTINUE 

EXTER=. FALSE. 

I  S=. FALSE. 

RAY=.  FALSE. 

CALL  PRF I1T( XD, SU ,XR , ATTENT , PHASET , ERV , ERVM , XFV , SFV ,.IN, RAY) 

ATT  E  N* ATTE  N +ATT E NT 

PHASE=PHASE+PHASET 

DB»ATTEN*8.6H58 

IF( PATH)  WRITE( 6 , 3U0)  ( XR( I ) , 1= 1 , 3) ,DB .PHASE 
IF(PATH)  RETURN 

IF( ERVM. LE. EPSA( I ) )  THEN 
PATH=.TRUE. 

END  IF 
C 

C  NEW  GUESS  FOR  INITIAL  DIRECTION 
C 

SS=0. 

DO  130  1=1,3 

SU(  I  )=DUM(  1  )+RC*ERV(  I  )/Sl, 

SS=SS+SU( 1 )**2 
130  CONTINUE 
SS=SORT(SS) 

DO  135  1=1,3 
SU(I)=SU{  D/SS 
135  CONTINUE 

I F ( PATH )  GO  TO  105 
IT=1T-H 

IF( IT.CT.NIT)  GO  TO  140 
GO  TO  105 
140  CONTINUE 

WR ITE( 6 , 150) 
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150  F0RMAT(1H  ,'IT=MAX,  NO  CONVERGENCE') 

300  F0RMAT(1H  ,5( 1PE10. 3, IX) ) 

RETURN 

END 

******** ************* ****** ***** *********** ******* ********** ********** * 
SUBROUTINE  PRFIEP 
C 

C  CREATES  FILE  OF  PLUME  INTERFERENCE  PROPERTIES  FORM  SPF 
C  FLOWFIELD  FILE  FF(I,J,K) 

C 

C  REVISED  FILE  HAS  THE  QUANTITIES  AS  FOLLOWS 
C  F( I , J , 1 )=RADIUS  AT  AXIAL  LOCATION  I,  RADIAL  LOCATION  J 
C  2  =INDEX  OF  REFRACTION 

C  3  “ATTENUATION  CONSTANT 

C  4  “PHASE  CONSTANT 

C  5  “WAVENUMBER* TURBULENCE  SCALE  LENGTH 

6  “SCATTERING  CROSS-SECTION 

7  “VELOCITY 

8  “TURBULENCE  INTENSITY 

9  “RADIAL  DERIVATIVE  OF  INDEX  OF  REFRACTION 
10  “AXIAL  DERIVATIVE  OF  INDEX  OF  REFRACTION 


INCLUDE  ' PRFICM.FOR/NOLIST' 

DIMENSION  GG( 10) ,HM(9) ,HP(9) 

DATA  RTHOM/ 2.81 33E- 1 5/ 

DO  50  1=1 ,NZ 
DO  40  J“1 , NPTS( I ) 

GG( 1 )=FF( I , J, 1 ) 

PFRE0=8. 97*SQRT( FF( I,J,4)) 

CALL  PRFIPC( BFREQ , PFREQ, FF( I , J ,6) ,GG( 3 ) ,GG( 4 ) ,GG( 2 ) ) 
GG(5)=2.*PI*GG(2)*BFREQ*FF( I,J,7)/C 
DENOM=( 1 .+( FF( I , J , 6)/ ( 2. * PI* BFREQ) )**2 )**2 
GG(6)“8.*( PI**3)* .0528*FF(1,J,5)*( ( RTHOM* FF( I,J,4))**2)* 
1  (FF( I , J, 7)**3)/DENOM 

GG(7)=FF( I, J,8) 

GG( 8)=FF( I , J , 9) 

DO  30  K=1 , 8 
FF(  I,J,K)=GG(K) 

30  CONTINUE 
40  CONTINUE 
50  CONTINUE 


COMPUTE  RADIAL  AND  AXIAL  DERIVATIVES  OF  INDEX  OF  REFRACTION  IN 
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DO  70  I»1,NZ 
DO  60  J-l,NPTS(l) 

ANR=0. 

BNR=0. 

DR>1“0. 

DRP“0. 

IF( J.GT.  1)  DRM-FK! £,J,1)-FF( I.J-l, 1) 

IF(  J.CT.  1  .AND.  DRM.NE.O.  )  BNR=(  FF(  i ,  J ,  "> )  -  F  F  (  I ,  J- 1 , 2  )  ) /DRM 
IF(J.LT.NPTS(I))  DRP=FF(I,J+1 , 1 )-FF( I , J , 1 ) 

IF( J . LT. NPTS( I ) . AND. DRP. NE. 0. )  ANR“! FF( I , J+ 1 , 2)-FF( I , J , 2 ) ) /DRP 

lF(J.EO.l)  BNR=ANR 

IF! J.EQ.NPTS!  I) )  ANR=BNR 

FF(  I ,  J,9)  =  . 6*( ANR+BNR) 

IF(L.GT.I)  CALL  PRFILK Z( l-l ) ,FF( I , J , l ) ,HM) 

IFU.LT.NZ)  CALL  PRFILK  Z(  1+)  ), FF(  I  ,J,  1 )  ,HP) 

I F ( L.CT. 1 )  BNZ=( FF( l , J , 2 )-HM( 1 ))/(Z( I)-Z( 1-1)) 

IF( I.LT.NZ)  ANZ-(HP( 1)-FF( I , J , 2) )/(Z( 1+1 )-Z( I ) ) 

IF(I.EQ. 1)  BNZ=ANZ 
IF(  L.EO.NZ)  ANZ-BN'Z 
FF( I , J , 10)=. 5*( ANZ+BNZ) 


60  CONTINL'E 
70  CONTINUE 
RETURN 
END 

*******************  V;  **************************** *********************** 


C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 


SUBROUTINE  PRF IFF 

READS  SPF  FLOWFIF.LD  FILE  INTO  ARRAY  FF(I,J,K) 


FF(l,J,l)  “  RADIUS! M)  AT  AXIAL  LOCATION  I,  RADIAL  LOCATION  J 

2  “  TEMPERATURE! K) 

3  =  PRESSURE! NT/ M** 2 ) 

It  =■  ELECTRON  NUMBER  DENSITY 

5  =  MEAN  SQUARE  ELECTRON  NUMBER  DENSITY 

6  =  CFREQ 

7  -  TURBULENCE  MACRO  SCALE 

8  “  MEAN  VELOCITY 

9  =  MEAN  SQUARE  TURBULENCE  VELOCITY 


INCLUDE  ' PRF ICM. FOR/ NOLI  ST' 
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DATA  NFF  /  7  /,  I  /  1  / ,  NZMAX.NRMAX  /  60,50  / 


OPEN(UNIT=NFF,  STATUS =’ OLD' ,  FILE=FILEFF,  READONLY) 

100  READ( NFF, 1000, END=300)  Z(I),  NPTS(I) 

D  WRITE(6 , 2001 )  I ,Z( I) ,NPTS( I) 

IF(NPTS( I) .GT.NRMAX)  THEN 

WRITE(6 ,* )  '  FLOWFIELD  EXCEEDS  RADIAL  DIMENSION.' 

STOP 

END  IF 

DO  200  J=1 ,NPTS( I) 

READ( NFF, 2000)  ( FF( I , J,K) ,K=1 , 9) 

D  WRITE(6 , 2002)  J,(FF( I ,J,K) ,K=1 ,9) 

200  CONTINUE 

1  =  1+1 

IF( I.LE.NZMAX)  GO  TO  100 
300  NZ  =  I  -  1 


1000  FORMAT(E 10. 3,110) 

2000  FORMAT ( 9E 1 0 . 3 ) 

2001  FORMAT( I5,E10.3,I5) 

2002  FORMAT(I5,9E10.3) 

RETURN 

END 

******************  ********  *******  ****  *****  *  ****  ***************  **  ******* 
SUBROUTINE  PRFIIN 

READS  IN  INITIAL  DATA  FOR  PRFIC 


INCLUDE  'PRFICM. FOR/NOLIST’ 
DATA  PRAD/1.745329252E-02/ 


C  CARD  1 

READ( 5,6)  I,  FILEFF 
FILEFF(I+1)=*0 

C  CARD  2 

READ( 5, 7)  TITLE 
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C  CARD  3 

READ( 5,1) 

C  CARD  4 

READ( 5,1) 

C  CARD  5 

READ( 5,1) 

C  CARD  6 

READ( 5,1) 

C  CARD  7 

READ( 5,1) 

C  CARD  8 

READ( 5,1) 

C  CARD  9 

READ( 5,1) 

C  CARD  10 

READ( 5,1) 

C  CARD  1 1 

READ( 5,1) 

C  CARL  12 

READ( 5,1) 

C  CARD  1 3 

•  READ( 5,1) 

C  CARD  14 

READ( 5,1) 

C  CARD  i 5 

READ( 5,1) 

C  CARD  16 

READ( 5,2) 

C  CARD  1 7 

READ( 5,2) 

C  CARD  18 


BFRE(),ZLL,ZUL 

(XT(J),J-1.3) 
(DCT(J),J-I,3) 
(DCT0( J) ,J“1 , 3) 
(VT(J),J-1,3) 
(XR( J) ,J=1 ,3) 
(DCR(J) ,J»1 ,3) 

(DCRO(J) , J-l ,3) 
(VR(J),J-1,3) 
(XM(J),J-1,3) 
(DCM(J) ,J=1 ,3) 
(DCMO(J) ,J=1 ,3) 
(VM(J),J-1 ,3) 

NCUTS , RC , NIT 

KPSA 
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READ( 5,5)  RAYTK , FSCAT 


READ  ANTENNA  GAIN  DISTRIBUTIONS 


IF( FSCAT)  THEN 
READ( 5,60)  NANG 
60  FORMAT! 2 J 5) 

IF(  NANG(  D.LE.O)  GO  TO  80 

READ( 5,70)  ( ANGLE ( 1 ,J) ,J=1 ,NANG( 1 ) ) 

READ( 5 ,70)  (GA IN( 1 ,J) , J= 1 , NANC( 1 ) ) 

80  CONTINUE 

IF( NANG( 2 ) . LE. 0)  GO  TO  90 

READ( 5,70)  ( ANGLE ( 2 , J ) , J= 1 , NANG( 2 ) ) 

READ( 5,70)  (GAIN(2,J) , J= I ,NANG( 2 ) ) 

GO  TO  130 
90  CONTINUE 
70  FORMAT! 10F8.0) 

IF(NANG( 1 ) . EQ.O. AND. NANG( 2) .EO. 0)  GO  TO  130 
IF(NANG(2).GT.O)  GO  TO  110 
DO  100  1=1 ,NANG( 1 ) 

ANGLE( 2 , I)=ANGLE( 1,1) 

CAIN(2,I)=CAIN(1,1) 

100  CONTINUE 

NANC(2)=NANG( 1) 

GO  TO  130 
110  CONTINUE 

DO  120  1= 1 , NANG( 2 ) 

ANGLE! 1 ,I)=ANGLE(2,I) 

GAIN( 1 ,I)=GAIN(2,I) 

120  CONTINUE 

NANG( 1)=NANG(2) 

END  IF 
130  CONTINUE 

1  FORMAT! 3F 10.0) 

2  FORMAT! 2F 10. 0,13) 

5  FORMAT! 2L1) 

6  FORMAT! Q.20A1) 

7  FORMAT! 18AA) 

DO  3  1=1,3 

DCT( I )=COS( DCT( I )*PRAD) 

DCT0( I )=COS( DCTO! I )*PRAD) 

OCR! I)=COS!DCR! I)*PRAD) 

DCRO( I )=COS( DCRO! I )*PRAD) 

DCM! I ) =COS ! DCM( I ) *  PRAD ) 
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DCMO( 1 )=COS(DCMO( I)*PRAD) 
5  CONTINUE 

C  OUTPUT  INITIAL  CONDITIONS 

CALL  DUMP  IN 

RETURN 

END 


SUBROUTINE  PRKIIC 
C 

C  READS  INPUT  CHANCES  FOR  MULTIPLE  CASE  RUNS  - 

C  ALLOWS  ONE  TO  INPUT  NEW  VALUES  FOR  VARIABLES 

C  ON  CARDS  3-1H. 

C 

C  INPUT  CHANCES  PROCEDURE  : 

C  l)A  CARD  CONTAINING  THE  TOTAL  NUMBER  Or  VARIABLES 

0  TO  BE  CHANGED  FOR  THE  NEXT  RUN. 

C  2)  CARD(S)  CONTAINING  THE  SEOUKNCE  OF  CHANGES  DESIRED. 

C  THE  FORMAT  OF  THIS  CARD(S)  IS  : 

CARD  NO.  OF  VARIABLE  TO  CHANGE .VARIABLE  NO.  ON  THIS  CARD, 
NEW  VALUE  FOR  THIS  VARIABLE  ... 

C  REPEAT  THIS  SEOUENCE  FOR  AS  MANY  CHANGES  SPECIFIED  ON 

C  CARD  1).  CARD  2)  MAY  BE  REPEATED  AS  MANY  TIMES  AS  NECESSARY 

C  TO  COMPLETE  THE  NUMBER  OF  CHANGES  SPECIFIED. 

C  FOR  EXAMPLE,  TO  CHANGE  THE  VARIABLE  ".7.UL"  ON  CARD  3  TO  THE 

C  QUANTITY  5000  AND  TO  CHANGE  THE  VARIABLE  " RAYTR"  ON  CARD  IB 
C  TO  THE  VALUE  FALSE,  THE  FOLLOWING  CARDS  SHOULD  HE  INPUT  ... 

C  2 

C  3,3, 5000 ,  l 8 , 1 , F 


INCLUDE  ' PRFICM. FOR/NOLIST' 

LOGICAL  DINP( 2) 

DIMENSION  CVC3.3:  IfO  ,  IC(6  5)  ,  IV(f>5) 
EQUIVALENCE  < BFREo,CV( 1.5)) 


116 


non 


NVJC  TP  6386 


EQUIVALENCE  ( RAYTR , DINP( 1 ) ) 
PRAD=1 . 745329252E-02 


READ( 5 ,* ,END=400)  NC 

1UO  READ( 5,* ,END=200)  ( IC(N) , IV(N) ,CV( IV(N) , IC( N) ) ,N=1 ,NC) 
200  DO  300  N= 1 , NC 

IF(IC(N).EQ.5  .OR.  IC(N).EQ.6  .OR.  IC(N).EQ.9 
<f  .OR.  IC(  N) .  EQ.  10  .OR.  IC(N).EQ.  13  .OR.  IC(N).EO.IA)  THEN 
CV( 1V( N) , I C ( N ) )=COS( CV( IV(N) ,IC(N))*PRAD) 

END  IF 

300  CONTINUE 

C  OUTPUT  INITIAL  CONDITIONS 

WRITEC6,*)  *1' 

CALL  DUMPIN' 

RETURN 

490  STOP 
END 


SUBROUTINE  DUMPIN 

DUMPS  INITIAL  DATA  TO  PRINTER 


INCLUDE  'PRFICM. FOR/NOLIST' 

BYTE  IDATE(9) 

DIMENSION  ACT( 3) ,ACR(3) ,ACM(3) ,ACT0(3) ,ACR0(3) ,ACM0(3) 
PRAD=1 . 74 53292 52E-02 
CALL  DATE ([DATE) 

DO  10  1=1,3 

ACT( I )=ACOS( DCT( I ) ) / PRAD 
ACR( I )=ACOS( DCR( I ) ) / PRAD 
ACM( I ) = ACOS ( DCM( I ) ) / PRAD 
ACT0( I )=ACOS( DCT0( I ) ) /PRAD 
ACR0( I )=ACOS( DCR0( I ) ) / PRAD 
ACM0( I )=ACOS( DCM0( I ) )/PRAD 
10  CONTINUE 
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WKITE(b,  100(1)  1DATE,  TITLE 
UK  I  I'E(b  ,  2()Oi  0  BFRK0/1.E9,  FILEFF 
UK ITE( h , 3000 )  XT,  <R ,  XM 
WRL TE( 6,4(11)0)  ALT,  ACR,  ACM 
WRlTE(b.5UiM>)  ACTO,  ACRU,  ACMO 
UK ITK(b , boon)  VT,  VR ,  VM 
WR I  I'E( 6 , 7000 )  RAYTR,  FSCAT 
wRITE(6,80nO)  NCl'TS.  RC,  EPSA,  NIT 

C  FORMATS : 

1000  F0RMATUH1,"/  50X,  ’  NAVAL  WEAPONS  CENTER,  Chi  NA  LAKE’// 

••■•39X, ’EXHAUST  PLUME  RADAR  FKEOUENCY  INTERFERENCE  CODE(PRFIC)’ 
■‘//38X, ’WRITTEN  BY  AERONAUTICAL  RESEARCH  ASSOCIATES  OF  PRINCETON’ 
••••//t>2X,9Al///45X,  18A4,  //) 

2000  FORMAT! 4 5X. 'TRANSMITTER  FkEolENCY  «  ’.IPEl'M,'  GHZ’/ 

•" 4 5X ,  ’  FLOW F  l ELD  FILE  -  ’  ,2oAl  ///) 

JOOO  FORMAT! 59X , ’ GEOMETRIC  DATA ’ / / 

«J  IX,  ’TRANSMITTER’  , 23X,  ’RECEIVER’  ,2hX, 

■•••’MISSLK’//  ’  POSITION,  M’  ,4X,  3! 4\,  9X,  * X’  , 9X ,  ’  Y  '  ,9X,  ’ Z  ’  ;/  •  bX, 

3(4X, 3(  1  PE10.  )))/) 

4 OUO  FORMAT! ’  DIRECTION  OF’ , 3X, 3! -.X, 5X, ’ALPHA’ ,  bX, ’ BETA’ , 3X, ’GAMMA’  )/ 
NORMAL,  DEC’,4X,3(4X,3!lPF.lo,3))/) 

5000  FORMAT!  ’  REFERENCE*  ,hX,3(4X,  5X,  ’ALPHA’  ,bX,  ’  BETA’  ,  3X,  ’GAMMA’  )/ 

'  DIRECTION,  DEC’ ,!X,3(4X,3(  IPEIo.3))/) 

6000  FORMAT!  ’  VELOCITY  ’,  7X,  3(  4X,9X,  ’  U  ’  ,9X,’V’  ,9X,’W)' 

«  ’  COMPONENTS,  M/SEC’  ,2X,3(  1PE10.  3)  ,  2!  4X,  3(  1  PEIO.  3)  ) .//  > 

7000  FORMAT! 25X, ’ RAYTRACE  CALCULATION  =  ’ , Ll/4  3X, ’ SCATTERING  ’, 
•■••’CALCULATION  -  ’,Ll/> 

8000  FORMAT! 2 5X, ’NO.  OF  STEPS  IN  INTEGRATION  THROUGH  PLUME  =  ’, 
l  I  PE  I".  )/ 

-.3X,  ’RELAXATION  CONSTANT  FOR  RAYTRACE  ITERATION  =  ’  ,  I  Ft  1 0.  3/ 

:f-*5X,  ’  ERROR  CRITERIA  FOR  RAYTRACE  CONVERGENCE  •  ’,1PE10.3,’,  ’, 

PE  in.  3/45X,’  NUMBER  OF  ITERATIONS  IN  RAY  TRACE  =  ’,13) 

RETURN 

END 

ft*************************************************************'*'*^****** 

SUBROUTINE  PRF I IT( XO , SO, XF , ATTEN, PHASE , ERV , ERVM , XFV , SFV , I N , RAY ) 

C 

C  COMPUTES  THE  PATH  FROM  THE  INITIAL  POINT  XO  TO  THE  FINAL  POINT 
C  XF  WITH  INITIAL  DIRECTION  SO 

C 

C  ATTEN- INTEGRATED  ATTENUATION  DUE  TO  ABSORPTION  AND  SCATTERING 
C  PHASE- INTEGRATED  PHASE  CHANGE 
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C  ERV=ERROR  VECTOR  BETWEEN  ACTUAL  END  POINT  AND  COMPUTED  END  POINT 
C  ERVM=MACN1TUDE  OF  ERROR 

C  XFV-COMPUTED  END  POINT 

C  SFV=COMPUTED  DIRECTION  OF  RAY  AT  END  POINT 

C  IN=T  ;  XO  AND  XF  INSIDE  THE  PLUME 

C  F  ;  OUTSIDE  THE  PLUME 

C  RAY=T  ;  DO  RAY  TRACE  CALCULATION 

C  F  ;  DO  LOS  CALCULATION 

INCLUDE 'PRFICM.FOR/NOL 1ST' 

LOGICAL  IN, RAY 

DIMENSION  XO(3),SO( 3),XF( I) ,ERV(3) ,XFV(3) ,SFV(3) ,X(3) ,SU(3), 

1  HOLD( 9 ) , HNEW( 9 ) , CURV( 3 ) 

C 

C  FIRST  SEGMENT  DOES  THE  LOS  CALCULATION  FOR  POINTS  OUTSIDE  THE 

C  PLUME 

C 

ATTEN=0. 

PHASE=0. 

SO-0. 

ERVM=0. 

DO  10  1=1,3 

SQ=SCH-(XF(I)-X0(I))**2 
10  CONTINUE 
SL=SORT( SQ) 

IF(IN)  GO  TO  20 

DO  15  1=1,3 

XFV( I )=X0( I )+SL*S0( I) 

ERV( I )=XF( I )-XFV( I ) 

ERVM=ERVM+ERV( I)**2 
SFV( I )=S0( I ) 

15  CONTINUE 

R=SQRT(XFV( 1 )**2+XFV( 2)**2) 

CALL  PRFILI(XFV( 3) ,R,HNEW) 

ERVM=SQRT( ERVM) 

PHASE=2.*PI*BFREQ*SL/C 

IF( PATH. AND. EXTER)  WRITE(6,200)  (XFV( I ) , 1=1 , 3) ,ATTEN, PHASE, 

1  (HNEW( I) , 1=1 , 7) 

RETURN 


SECOND  SEGMENT  DOES  LOS  CALCULATION  WITHIN  THE  PLUME 


20  CONTINUE 


119 


NWC  TP.  6336 


D'-SL/ NCUTS 
h: i»0. 5*DS 
D  JO  1-1,3 
>:  :  >»xu(  i ) 

S.  L)=SO(I) 

30  C  NTI.NUE 

R».-ORT(X(  1  )**2+X(  2  )**2) 

CALL  PRFILK  X(  3)  ,R,IIOLD) 

i  F  »  i . 

I F  RAY)  GO  TO  70 
40  C  NTINUE 
DC  30  1-1,3 
■  X  :)-X(I)+SU(I)*l)S 
50  C  NTINUE 

:<«>ORT(X<  1  )**2+X(2)**2) 
call  prfiliixu)  ,r,hnew) 

A I7EN-ATTEN+. 3*( HOLD! 2 )+H0LD( 5 )*SCAT( UOLD( 4 ) )+HNEW( 2 )+HKEW( 5 ) * 
1  SCAT! HNEW( 4 ) ) )*DS 

F-ASE-PHASE+. 5* ( HOLD! 3 )+HNEW( 3 ) ) *DS 
D  60  J=l,9 
K  "_D(  J  ) *HNEW(  J  ) 

60  C  .  NTINUE 
S5  =  SS*-DS 

dif-abs(sl-ss) 

:i-ATTEN*8.6858 

C  WRITE  STATEMENT  FOR  DETAILED  DISTRIBUTION  ALONG  FINAL 

C  INTEGRTION  PAIH 

C  IF  PATH)  WRITE( 6, 200)  ( X( I ) , I- 1 , 3) , DB .PHASE, ( HNEW( I ) , I- 1 , 7 ) 

IF  SS. LT. SL. AND. DIF.GT. HDS )  GO  TO  40 

EEVM-O. 

r:  65  1-1,3 

yjv(i)-xd) 

ERV( I)-XF( I )-X( I) 

SFV(I)-SOU) 

65  C  :  NTINUE 

ERVM-SCALP!  ERV ,  ERV) 

ERVM-SQRT(ERVM) 

RETURN 

C 

C  THIS  SEGMENT  DOES  THE  CALCULATION  WITHIN  THE  PLUME  WITH  RAY 
C  CVRVATURE 
C 


70  C  NTINUE 
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ERVM=SL 
80  CONTINUE 
ERVMO=ERVM 

R=SQRT(X( 1 )**2+X( 2 )**2 ) 

C0ST=X( I ) / R 
SINT=X(2)/R 

TD0TGN=SU(  l)*HOLD(8) *C0ST+SU( 2 )*H0LD( 8) *S INT+SU( 3)*HOLD( 9 ) 

Ol'RVC  1  )  =  (  HOLD(  8)*C0ST-SU(  1  )*TDOTGN)/HOLD(  1 ) 

CURV( 2 ) =( H0LD( 8 )*S INT-SU( 2)*TDOTGN)/HOLD( 1 ) 

CURV (3)=(UOLD(9)  -SU( 3)*TD0TGN)/H0LD( 1 ) 

DO  90  1=1,3 

X(  I  )=X( I )+SU( I )*DS+. 5*CURV( I ) *DS*DS 
SU( I)=SU( I)+CURV( I)*DS 
90  CONTINUE 

R=SORT(X(  1 )**2+X( 2 )**2) 

CALL  PRFILI(X( 3) , R,HNEW) 

ATTEN=ATTEN+.5*(HOLD(2)+HOLD( 3)*SCAT(HOLD(4) )+HNEW( 2)+t(NEW( 5 )* 

1  SCAT(HNEW(4)))*DS 

PHASE=PHASE+. 5*( H0LD( 3)+HNEW( 3) )*DS 
SS=SS+DS 
DO  100  J=1 , 9 
H0LD( J ) *HNEW ( J ) 

100  CONTINUE 
ERVM=0. 

DO  110  1=1,3 
ERV( I )=XF( I )-X( I ) 

ERVM=ERVM+ERV( I )**2 
110  CONTINUE 

ERVM=SQRT(ERVM) 

DB=ATTEN*8. 6858 

WRITE  STATEMENT  FOR  DETAILED  DISTRIBUTION  ALONG  FINAL 
INTEGRTION  PATH 

IF( PATH)  WRITE( 6 , 200)  (X( I) , T=1 , 3) ,DB, PHASE, (HNEW( I) , 1=1 , 7) 

IF( ERVM.LT. ERVMO. AND. ERVM.GT. HDS)  GO  TO  80 
DO  120  T=1 , 3 
XFV( I )=X( I ) 

SFV( I )=SU( I ) 

120  CONTINUE 

200  FORMATC1H  ,  1 2( 1PE10. 3 , IX) ) 

RETURN 

END 

***************** ********************************** *********** ******** 

SUBROUTINE  PRFILI( ZL, RL, VZR) 

GIVEN  AXIAL  LOCATION,  ZL,  AND  RADIAL  LOCATION,  RL, 
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C  PRFILI  LINEARLY  INTERPOLATES  TO  FIND  THE  FLOWF1ELD 

PROPERTIES  AT  ZL,  RL  AND  RETURNS  THEIR  VALUES  IN  VZR. 


;.NCLUDE  ’PRFICM.  FOR/ NOLIST' 

DIMENSION  VZR(9),  DATZ(2:10),  DATXP(2:10) 

DATA  VZRAMB  /  i.,  8*0.  /,  I LAST  111 

YZRAMB(  ))  =  2.*PI*BFRE0/C 
YZRA>1B(6)  =  FF(NZ,NPTS(NZ) ,7) 


C 

C  SEARCH  FOR  BOUNDING  AXIAL  STATIONS 

c 

!F( (ZL.LT. Z( I ))  .OR.  (ZL.GT.Z(NZ)))  GO  TO  600 
DO  100  I  * I  LAST , NZ 
loti  IF(ZL.LE.Z(  H)  GO  TO  20n 

200  1  =  1-1 

C  I LAST  =  I 

DO  400  IZ*1  ,2 
DO  300  J=2,NPTS( I) 

300  :f(abs(rl).lt.ffu,j,  i))  go  to  3 1 o 

CALL  SFVMV( VZRAMB , DATZ , 9 ) 

GO  TO  330 

310  RRAT  =  (ABS(RL)~FF( I , J- 1 , l))/(FF( I , .1 , 1 )-FF(  I ,  J- 1 , 1)) 

DO  >20  K=2 , 1 0 

320  DATZ(K)  =  FF(l.J-l.K)  +  RRAT*( FF( I , J,K)-FF( 1 ,J-1 , K) ) 

C 

C  SAVE  RADIALLY  INTERPOLATED  VALUES  A I  THE  LOWER  BOUNDING 

C  AXIAL  STATION  AND  REPEAT  FOR  UPPER  BOUNDING  STATION. 

C 

330  1F( IZ.EQ. 1 )  CALL  SFVMV ( DATZ , DATZP , 9 ) 

400  1=1+1 


C 

C 

C 


300 


600 


INTERPOLATE  BETWEEN  AXIAL  STATIONS 
1  =  1-1 

Z RAT  =  (ZL-Z( I— 1))/(Z(I )— Z( I — I ) ) 

DO  500  K=  2 , 10 

VZR(K-i)  =  DATZP(K)  +  ZKAT*( DATZ( K )-DATZP( K) ) 

RETURN 

CALL  SFVMV( VZRAMB, VZR  ,9) 

RETURN 

END 

'If+rk'k'k'k+ck'k-k’k'k-kit'k'k'k'k'k-k'k'k'k'k'k-iK’tz'k'k’k'k'k-k'k'k'k'k'k-ti'k'k'k'k-itifk'k'klfk’k'k'k'k'k'k’k  kii-k-kb-tyk 
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PROGRAM  PRFIMM 

MAIN  CALLING  PROGRAM  FOR  PRFIC 

INCLUDE  ' PRFICM. FOR/NOLIST' 

CALL  PRFIIN 
CALL  PRFIFF 
CALL  PRFIEP 
5  CONTINUE 

IF( ZLL. LE.O. )  ZLL=Z( 1 ) 

IF( ZUL. LE. 0. )  ZUL=Z( NZ) 

WR  ITF,(  6 , 20)  ZLL, ZUL 
WRITE(6,30) 

CALL  PRF l DT( ATTEN , PHASE , DOPP ) 

DPHASE=0, 

DO  15  1=1,3 

DPHASE=DPHASE+(XR(l)-XT(I))**2 
15  CONTINUE 

DPHASE=2 .  *P1*BFRF.Q*S0RT(  DPHASE)  /  C 
DPRASF.=PHASE-DPHASE 
ATTEN=ATTEN*8. 6858 
DOPP=DOPP/ 1000. 

WRITE(6, 10)  INSECT, ATTEN, PHASE, DPHASE, DOPP 

IF(FSCAT)  THEN 
WRITE(6, 100) 

IF(NANG( D.NE.O)  WRITE(6,110)  ( ANGLE ( 1 , J ) , J= 1 , NANG( 1 ) ) 
IF(NANGd).EQ.O)  WRITE(6,150) 

WRITE(6, 120) 

IF(NANG( D.NE.O)  WRITE(6,110)  (GAIN( 1 , J) , J=1 ,NANG(  1 ) ) 
WRITE(6, 130) 

IF( NANG( 2) .NE.O)  WRITE(6,110)  (ANGLE( 2 , J) , J= 1 ,NANG( 2) ) 
IF(NANG(2).EQ.O)  WRITE(6,150) 

WRITE(6, 140) 

IF( NANG (2). NE.O)  WRITE(6,110)  (GAIN( 2 , J) , J= 1 , NANG( 2 ) ) 
CALL  PRFISC 
END  IF 

CALL  PRF IIC 
GO  TO  5 

10  FORMAT!  1H  ,///  40X,  ’  INTERSECTION  WITH  PLUMF.=  ’  .L1/52X, 
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1  '  ATTENUATION-  '.FI0.4,*  DB  7 

2  --X,  'PHASE-  ', IPE12.5/47X, 'PHASE  DIFFERENCE-  '.1PE12.5/ 

3  T-<x, 'DOPPLER  SHIFTED  FREQUENCY-  ' .1PE10.3,'  KHZ') 

20  F  RMAT(  1H  ,44X,'PLi:ME  LENGTH  BETWEEN"  ,  1 PE10.  3 , '  AND  '  .LPE10.3, 

i  '  'VIII) 

30  F'.RMAT/  1H1 ,44X,  'PROPERTIES  ALONG  RAY  TRAJECTORY’// 

1  -X, ' X,M' ,8X, 'Y,M' ,HX, 'Z,M* , 3X, ' ALPHA, DB ' , 5X, ' BETA*L' , 1UX, *N' , 

2  2X, 'ALPHA, I/M' ,3X, 'BETA, I/M’ , 4X, ’ K*LAMDA’ ,2X, ' SIGMA, l/M' , 

3  V\,  'U,M/S'  ,3X,'(D/U)**27) 

100  F  oRMAT/  I H 1 ,  /  /  /  44X ,  'ANTENNA  GAIN  PATTERNS’  /  / //44X, 

1  ' TRANSMITTER’ ///IX, ' ANGLES , DEG : ' ) 

110  FORMAT/ l H  , 1 5X, lo( FIO. 2) / ( 1 6X, 10(  F 10. 2) ) ) 

120  FORMAT/  III  , /  /  3X,  'GAIN,  DB:’) 

l  50  r  ■■•■LMAT/IH  ,/ / / /  1  2X , '  RE CE I VER  ’  /  / /  l  X,  ’  ANGLES, DEG:  '  ) 

140  FORMAT/ 1 H  ,//3X,'GAIN,  DB:') 

130  F>RMAT(1H  ,20X,  'UNIFORM  DISTRIBUTION',  GAIN-0.  FOR  ALL  ANGLES’) 

END 

•k'k'k'k'k'k'k'kii'k'kir'k’k-k'k'k'k'k-k'k'k'k'k'k'k'k'k'k'k'k'k'k'k’k'k'k'k-k'M-k'k’k'k'kic'k'k^’kic'k'k'k'kit'M'k'k'k^it'k-k^cit'k-k'k'kit 

SUBROUTINE  PR  F I  PC  (  B  F  R  EO ,  PF  R  E<) ,  C  F  R  EO ,  A  LPHA ,  BET  A ,  R I ) 

C  CALCULATE  ATTENUATION  AND  PHASE  CONSTANT 

C  BFREQ-PROPAGAT ION  FREQUENCY ,  HZ 

C  PFREO-PLAsMA  FREOUENCY.HZ 

C  C r RE 0=F. LECTRON-NE UTRA L  COLLISION  FREQUENCY,  1/SEC 

ALPHA-ATTENUATION  CONSTANT, 1/M 
BETA  -PHASE  CONSTANT, 1/M 
R I- INDEX  OF  REFRACTION 

REAL  KR , K I , KMAC 

DATA  C, PI/ 2 . 998E+08 , 3. 141392653/ 

Rl-PFREQ/ BFREO 
WPROP-2 .  *  P  I*BFRE() 

R2-CFREQ/WPR0P 
R1S=R1*R1 
R2S=R2*R2 

KR-1.-R1S// 1.+R2S) 

IF/KR.LT.O. )  WRITE/6,  l) 

KI=RlS*R2/( 1 .+R2S) 

KMAG-SORT/ KR*KR+KI*K1 ) 

ALPHA-SORT/  .  5*  (  KMAC-KR)  )  *WPR()P/ C 
RI=SORT( . S*(KMAG+KR) ) 

5ETA=RI*WPR0P/C 
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RETURN 

1  F0RMAT(1H  PLUME  IS  LOCALLf  OVERDENSE’) 

END 

** ****** *★★****★* ***★*★*★ ************************ ************** ****** ** 
S UBROUT INE  PRF IP I ( XO , XF , XI , XE , MTYPE , NTYPE ) 

COMPUTE  THE  POINT(S)  AT  WHICH  THE  LOS  FROM  XO  TO  XF  INTERSECTS 
THE  PLUME  BOUNDARIES 


C  MTYPE. LT.O  :  XF-INITIAL  DIRECTION  COSINES 
C  MTYPE. GE.O  :  XF=COORDINATES  OF  FINAL  POINT 

C  NTYPE= 1 ,  LOS  PASSES  COMPLETELY  THROUGH  THE  PLUME  (TWO  SOLUTIONS, 

C  XI, AND  XE,  OR  NO  INTERSECTION,  INSECT=FALSE) 

C  NTYPE=2 ,  XO  OUTSIDE  THE  PLUME,  XF  INSIDE  PLUME  (ONE  SOLUTION,  XI) 

C  NTYPE=3,  XO  INSIDE  THE  PLUME,  XF  OUTSIDE  PLUME  (ONE  SOLUTION,  XE) 

C 

INCLUDE  'PRFICM. FOR/NOLIST’ 

LOGICAL  ROOT1 

DIMENSION  X0( 3) ,XF(3),XI(3) ,XE(3) ,XD(3) ,SU(3) ,ROOT(2) 

INSECT*. FALSE. 

ROOT 1=. FALSE. 

IF( MTYPE. GE.O)  GO  TO  5 
DO  1  1*1,3 
SU( I)*XF(  I) 

1  CONTINUE 
GO  TO  20 
5  CONTINUE 
SS*0. 

DO  10  1=1,3 
XD( I)=XF( I)-XO( I) 

SS=SS+XD( I)**2 
10  CONTINUE 
SL=SQRT(SS) 

DO  20  1-1,3 
SU( I)=XD( I)/SL 
20  CONTINUE 
NZM1=NZ-1 
DO  90  1=1 ,NZM1 
IP=I+1 

SLOPE=(FF( IP,NPTS( IP) , 1 )-FF( I ,NPTS( I),1))/(Z( IP)-Z( I) ) 

A=FF( I ,NPTS( I) , 1 )+SLOPE*(XO( 3)-Z( I) ) 

B=SLOPE*SU(  3) 

XP=X0( 1 )*SU( 1 )+X0( 2 )*SU( 2 ) 

XP5=X0( l)**2+XO(2)**2 
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ops=SU(  1  )**2+SU(2)**2 

AR=B* B-CPS 

.5R=A*B-XP 

CP.=A*A-XPS 

l:sc=br*br-ar*cr 

::-'(DlSC.LT.O.)  GO  TO  90 

IHAR.EO.il.)  GO  TO  60 

d:sc=si}rt(disg) 

R  t  =  ( -BR-DT  SC) / AR 
R2  =  ( -BR+D ISC ) / AR 

1F(RI.LT.0. .AN0.R2.LT. 0. )  GO  TO  90 
INSECT* . TRUE . 

ROOT( 1 ) =AM I N 1 ( R 1 , R2) 

IHRl.l.T.o. )  R00i'(  1  )=R2 
I K(  R2.LT.0. )  Ro(VP(  l  )=R  | 

:'.HX0(  ))+ROOT(  i)*SU(  3> 

1  F(  Z  1 . 1,T.  Z(  1 ) .OR.Z 1 .GT. Z(  IP))  GO  TO  40 
DO  30  J=  1 , 3 

XI  ( J  )  =  X0(  J  )+ROOT(  1  )*SUU> 

30  CONTINUE 
RoOTl  =.TRUK. 

IF (.  NTYPE.  El).  2  )  RETURN 
40  CONTINUE 

ROOT( 2 )=AMAX  U  R I ,R2) 

I F ( ROOT I  2 ) . NK. ROOT( 1 ) )  GO  TO  44 
ROOT(2)=ABS(Z(NZ)-XO( 3)) 

Go  TO  45 
44  CONTINUE 

Z>XU(3)+ROOT(2)*SU(3) 

IF U-2.LT. Z( I).0R.Z2.GT.Z( IP))  GO  TO  90 
4  j  CONTINUE 
1)0  50  J  -  l  ,  3 

XE(  J  )=XO(.l)+ROOT(  2  )*Sl!(  J ) 

50  CONTINUE 

IF(  NTYPE.  El).  3)  RETURN 
IF( ROOT l )  RETURN 
GO  TO  90 
60  CONTINUE 
P.  I  *  ~  1 . 

I F ( BR . NE . 0 , )  R1=-.5*CR/BR 
I r ( R 1 . LT.O. )  GO  TO  90 
INSECT*. TRUE. 

I F( NTYPE. EO. 3)  GO  TO  70 
RO0T( 1 )*R1 
DO  65  .1=1,3 

XI(  J)=XO(.t)+ROOT(  1  )*SU(J) 

r> 5  continue 
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ROOT >. TRUE. 
IF(NTYPE.EQ.2)  RETURN 
GO  TO  90 
70  CONTINUE 
R00T(2)=R1 
DO  75  J=1 ,3 

XE(J)=XO(J)+ROOT(2)*SU(J) 
75  CONTINUE 

IF(NTYPE.EQ.3)  RETURN 
IF(ROOTl)  RETURN 
90  CONTINUE 


IF( .NOT. INSECT)  RETURN 
IF( NTYPE. EQ. 3)  GO  TO  105 

IF(Z1.GE.Z(NZ).AND.SU(3).NE.0.)  ROOT( 1 )=(Z(NZ)-XO( 3) )/SU( 3) 

DO  100  1*1,3 

XI ( I )=X0( I )+R00T( 1 )*SU( I ) 

100  CONTINUE 

IF(NTYPE.EQ.2)  RETURN 
105  CONTINUE 

IF(Z2.GE.Z(NZ) .AND. SU( 3) .NE.O. )  ROOT( 2)=(Z(NZ)-X0( 3) )/SU( 3) 

DO  110  >1,3 

XE( I)=XO( I)+ROOT( 2)*SU( I) 

110  CONTINUE 
RETURN 
END 

*★*★*★**★★****★★**** ********************************************* ****** 


SUBROUTINE  PLOT(NO,A,N,M) 

INTEGER* 2  YLAB( 50) 

DIMENSION  OUT( 101 ) ,YPR( 1 1) ,ANG(9) ,A( 1 ) 

DATA  MX, MY/ 6, 5/ 

DATA  YLAB/13*'  ','F  ’,'R  ','E  ',’Q  ','U  \'E  ’,’N  ','C  ','Y  ’, 

It’  ' , '  D  ','1  ','F  ','F  '  , '  E  ','R  '  , '  E  '  , '  N  ','C  ','E,',’  ’, 

#’K  ' ,'H  ' ,’Z  ’  ,13*'  '/ 

DATA  ANG( 1 )/'*'/ ,  BLANK/'  '/ 

1  FORMAT  ( ' 1 ' , 60X, 7H  CHART  ,13,//) 

2  FORMAT  (IX, A2, IX, FI  1 . 1,5X, 1 0 1 A 1 ) 

3  FORMAT  ( 1X,A2) 

5  FORMAT(lOAl) 

7  FORMAT (  16X, 1 0 1 H .  .... 

1  .  .  .  .  .  .) 

8  FORMAT  (//.9X.11F10.  1) 

9  FORMAT  (// ,53X, '  SCATTERED  AMPLITUDE, DB/ HZ  ') 


NLL=50 

C  PRINT  TITLE 

20  WRITE (MX, 1 )NO 
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C  FIN’D  SCALE  FOR  BASK  VARIABLE 

XSCAL=( A( N)-A( 1 ) ) / ( FLOAT( NLL- 1 ) ) 

C  FIND  SCALE  FOR  CROSS-VARIABLES 

Ml=N+l 
M2=M*N 
YMAX=A(M1 ) 

DO  40  J=MI ,  M2 
IF(A(J)-YMAX)  40,40,30 
30  YMAX=A( J ) 

40  CONTINUE 

YMIN=YMAX-40. 

YSCAL= ( YMAX-YM I N ) / 1 00. 0 
C  FIND  BASE  VARIABLE  PRINT  POSITION 

XB=A( 1  ) 

L=1 

MYX  =  M-l 
1  =  1 

45  F=l-1 

XPR=XB+F*XSCAL 
[F( A(L)-XPR)  50,50,70 
C  FIND  CROSS-VARIABLES 

50  DO  55  IX* 1,101 
55  OUT( IX)=BLANK 
DO  60  J-l.MYX 
LL=L+J*N 

I F ( A( LL ) . LT . YM I N )  GOTO  60 
JP=((A(LL)-YMIN)/YSCAL)+l.O 
OUT! JP)=ANG( J ) 

60  CONTINUE 

C  PRINT  LINE  AND  CLEAR,  OR  SKIP 

WRITE(I  IX, 2 )YLAB( I ) , XPR , ( OUT( IZ) , 12=1,101) 

L=L+l 
GO  TO  80 

70  WRITE(MX, 3)  YLAB( I ) 

80  1=1+1 

IF( I -NLL) 4 5, 84, 86 
84  XPR=A( N) 

GO  TO  50 

C  PRINT  CROSS-VARIABLES  NUMBERS 

86  WRITE! MX, 7) 

YPR! 1 )=YM I N 
DO  90  KN=1 ,9 

90  YPR(KN+1 )=YPR(KN)+YSCP 4*10.0 
YPR( 1 1 )=YMAX 

WRITE! MX, 8) (YPR( IP) , IP=1 , 1 1) 

WRITE (MX, 9) 

RETURN 

END 

A***********************  ******  A*  ***:':*********  *******  */r****/r*******^-**^^r 
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SUBROUTINE  PRFIRZ( ZP , RP) 

C 

C  COMPUTES  PLUME  BOUNDARY  RADIUS  R  AT  AXIAL  LOCATION  2 
C 

INCLUDE  ' PRFICM.FOR/NOLIST* 

IF(ZP.LT.ZLL.OR.ZP.GT.ZUL)  GO  TO  30 
1=1 

10  CONTINUE 

IF(ZP.LE.Z( 1+1) .AND. ZP.GE.Z( I))  GO  TO  20 
1=1+1 

IF( I.GT.NZ)  GO  TO  30 
GO  TO  10 
20  CONTINUE 

RP=FF( I,NPTS( I),  1)+ 

1  (ZP-Z( I))*(FF( l+l ,NPTS( 1+1 ) , 1 ) — FF( I,NPTS(  I) , 1 )  )/(Z( [+1 )-Z(  I)  ) 
RETURN 
30  CONTINUE 

WRITE(6,40)  ZP, RP, ZLL.ZUL 

40  FORMAT( 1H  , ' NO  BOUNDARY  SOLUTION 2X, 4( 1  PE  1 0. 3 ) ) 

RETURN 

END 

a********************************************************************** 

SUBROUTINE  PRFISC 
C 

C  COMPUTES  VOLUME  SCATTERING  CONTRIBUTION  TO  ENERGY  RECEIVED  AT 

C  THE  ANTENNA 

C 

C  USES  GAUSSIAN  QUADRATURES  FOR  INTEGRATIONS  OVER  CIRCUMFERENCE, 

C  RADIUS,  AND  AXIAL  LENGTH 

C 

C 


INCLUDE  ' PRFICM. FOR/ NOLI ST' 

LOGICAL  IN, RAY 

DIMENSION  RTHETA( 10) ,WTHETA( 10) ,RRAD( 10) ,WRAD( 10) ,RZ(24) ,WZ(24) , 

1  X( 3) , ST( 3) , SR( 3) ,H( 9) ,DUM( 3),XI(3),XE(3), ERV( 3) ,XFV( 3) , SFV( 3) , 

2  SN(3) ,VP(3) ,SD(3) 

DIMENSION  VOLTH( 500) , VOLR( 500) ,VOLZ( 500) ,FVOLZ( 102) 

DATA  TWOPI/6. 283185307/ 

DATA  NNZ,NNR,NNTH,NF/24, 10, 10,50/ 

DATA  RTHETA/. 97390653, ,86 506337,. 67940956, .43339539, .14887434, 
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1  1 4887434, -.43 !  I9S  59,-.<,/9.4o9S6,-.86 506337, -.97390653/ 

DATA  UTHETA/ .06667  I  14  , .  1  494 5 1  35, .  2  1 908636  , .  26926672 , .  29532422  , 
1  . 2 9352 4 2 2 , . 269266  72 ,.2  1 9086  36 , .  1 494  5 13  5 , . 06667 1  34' 

DATA  RRAD/. 07652652, .22778585, .37370609, . 5 1 086700 , . 63605368 , 

1  .74633191 ,.8391 1 697 , .9 1 223443 , .96397 193 ,. 993 12860/ 

DATA  WKAD/ .  152753  39, . 1 49 1  7299 , . 1 4  2096 1 1 , . 13168864,.  11819453, 

1  .1019 1012, .08327674, .06267205, .04060 1 43 , .0 1 76 1 401 / 

DATA  RZ/ .001229  , .0064699, .0158755, .02941 85, .04/01 2  1, .0686 1  14, 

1  .0941209, .123428  ,. 1 5o4 1 88 ,. 1 9293  18, .23284 1  ,.2759659, 

2  .3221277, .  171 1 326, .4227753, .476834, .533097  ,.5913136, 

3  .6512442, .7126376, .7752363, .8387777, .9029953, .9676(99/ 
DATA  WZ/.0031 53 35, .00732755, .01 147723, .01 55793,. 01961616, 

1  .023570,6, .02742651 , .031 1 672 3, .03477722, .03824 1 35, 

2  .04154509, .0446745:), .04761666, .05035904, .05289019, 

3  .05519950, .0572/729, .0591 1 484 , ,o607o444, 

4  .0620 3942  ,.  063 1  1 4  1 9 ,. 0639242  ,  .06 44*>6 16,  .064 7 3770/ 


PATH=. FALSE. 

DOPPLER  SHIFT  PARAMETERS 

VTM=SCALP( VT , VT ) 

VRM»SCALP(VR,VR) 

VPM=FF( 1,1,7 )+VM( 3 ) 

VTM=SORT( VTM) 

VRM-St)RT(VRM) 

C  ESTIMATE  OF  MAXIMUM  FREQUENCY  SHIFT 

D  D  F = B  F REQ* ( VTM+V  RM+  2 . *  V  PM ) / C 

DFO=-DDF 

DNF=NF 

DELDF=2 . *DDF/DNF 
NFP=NF+1 

C  OUTER,  Z-INTECRAL 

C 

ALAM=C/ BFREQ 

CONST=(ALAM**2)/((2.*TWOPI)**3) 

ZNORM=ZUL-ZLL 

AEXP= 1 1 . / 6 . 

V0LZZ=0. 

DO  600  L= 1 , NFP 
VOLZ(L)=0. 

600  CONTINUE 
C 

DO  300  1=1, NNZ 
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MIDDLE,  R-INTEGRAL 

V0LZR=0. 

DO  610  L=1 ,NFP 
V0LR(L)=0. 

610  CONTINUE 

DO  200  J= 1 , NNR 
R=RRAD(J)*RPZ 

INNER,  THETA- INTEGRAL 

VOLZA=0. 

DO  620  L= 1 , NFP 
VOLTH(L)=0. 

620  CONTINUE 

CALL  PRFILI ( X( 3 ) , R, H ) 
IF(H(5).LE.0.)  GO  TO  100 
DO  100  K=  1 ,  liNTH 
THETA=PI*(  l.-RTHETA(K)) 

X( 1)=R*C0S(THETA) 
X(2)=R*SIN(THETA) 

CHECK  FOR  ADMISSIBLE  DIRECTIONS 

SST=0. 

SSR=0, 

DO  50  L-1,3 

ST(L)=X(L)-XT(L) 

SR(L)=XR(L)-X(L) 

50  CONTINUE 

SST=SCALP( ST, ST) 

SSR=SCALP( SR , SR) 

SST=SQRT( SST) 

SSR=SQRT( SSR) 

DO  60  L=1 , 3 
ST(L)=ST(L)/SST 
SR(L)=SR(L)/SSR 
60  CONTINUE 

COSPT=SCALP( DCT, ST) 
COSPR=-SCALP(  DCR, SR) 
IF(COSPT.LE.O. )  GO  TO  100 
IF( COSPR. LE. 0. )  GO  TO  100 
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ATTENUATION  AND  PHASE  KOK  THK  INCIDENT  AND  SCATTERED  KAY 
C 

C  INCIDENT  RAY 
C 

IT=  1 

70  CONTINUE 

CALL  PRK I P I  ( XT ,  ST,  X I , DOM ,-1,2) 

IN-. FALSE. 

RAY=. FALSE. 

CALL  PRF I IT( XT, ST, XI , ATTEN I , PHASE  1 , ERV . ERVM ,XFV , SFV , I S, RAY ) 
I  N= .  TRUE . 

RAY=RAY1R 

CALL  PRF II T( X I , ST , X , ATTEN? , PHASE 3 , ERV , E RVM , XEV , SFV , I N , RAY ) 
IF(RAYTR)  THEN 

I  F(  ERVM.  LE.  F.PSA(  2  ) )  CIO  TO  86 
SS=0. 

DO  75  I.»l,3 

XI  ( L)-XI( L)+RC*ERV( L) 

ST(L)=XI( L)-XT( L) 

SS=SS+ST(L)**2 
75  CONTINUE 
SS-SORT(SS) 

DO  80  L-1,3 
ST(L)=ST( L) / SS 
80  CONTINUE 

C  ATTENUATION  SET  TO  URGE  VALUE  IF  NO  CONVERGENCE 
IT-IT+1 

IF(  IT. (IT. NIT)  ATTF.N2=40. 

IF(  IT.GT.NIT)  GO  TO  85 
CO  TO  70 
END  IF 
85  CONTINUE 
C 

C  TRANSMITTER  GAIN 
C 

COSTA-SCALPI ST, DCT) 

IF(  COSTA.  LT.  - 1 .  )  COSTA— 1. 

IF( COSTA. GT.  1.)  COSTA-  1. 

TAANG-ACOS (  COSTA) *  1  80 .  /  P I 
EXPTG-0. 

IF(TAANG.GE.90. )  GO  TO  87 
CALL  PRFIAGCTAANG.TAG, 1) 

EXPTG- 10.**(TAG/10.) 

C 
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87  CONTINUE 

ATTENT-ATTEN1+ATTEN2 
PHASET-PRASE 1+PHASE2 
C 

C  SCATTERED  RAY 
C 

IT-1 

CALL  PRF I P I  ( X ,  SR ,  DIJM ,  XE  ,  - 1 , 3 ) 

90  CONTINUE 
IN-. TRUE. 

RAY-RAYTR 

CALL  PRF I IT( X , SR , XE , ATTEN 1 , PHASE 1 , ERV , ERVM , XFV , SFV, IN, RAY) 
IN=. FALSE. 

RAY=. FALSE. 

DO  95  L-1,3 
SR(L)=SFV(L) 

XE( L)=XFV( L) 

95  CONTINUE 

CALL  PRF I IT( XE , SR , XR , ATTEN2 , PHASE2 , ERV , ERVM , XFV , SFV , IN , RAY ) 
IF(RAYTR)  THEN 

IF( ERVM. LE. EPSA( 1 ) )  GO  TO  105 
SS-O. 

DO  110  L-1,3 
XE(L)=XE(L)+RC*ERV(L) 

SR(L)-XE(L)-X(L) 

SS=SS+SR(L)**2 
110  CONTINUE 
SS-SQRT(SS) 

DO  115  L-1,3 
SR' L)=SR(L)/SS 
115  CONTINUE 

C  ATTENUATION  SET  TO  LARGE  VALUE  IF  NO  CONVERGENCE 
IT-IT+1 

IF( IT.GT.NIT)  ATTEN2-40. 

IF( IT.GT.NIT)  GO  TO  105 
GO  TO  90 
END  IF 
105  CONTINUE 
C 

C  RECEIVER  GAIN 
C 

COSRA— SCALP(  SFV ,  DCR) 

IF(C0SRA.LT.-1.)  COSRA— 1. 

IF( COSRA. GT.  1.)  COSRA-  1. 


133 


o  n 


NWC  TP  6386 


RAANC=ACOS( COSRA)* 1 80. / P I 
KXPRG=0. 

1F(  RAANG. OK.  90.  )  CO  TO  107 
CAi.L  PRF I AG(  RAANG ,  KAC ,  2 ) 

EXPRG= 1 0.**( RAG/ 10.) 

C 

107  continue 

A TTENR=ATTKN 1 +ATTEN2 
PHASER=PHASE 1+PI1ASE2 
PHASE-PHASKT+PHASER 
ATTEN-ATTENT+ATTENR 
COSSC=SCALP( ST,  SR) 

I F ( COSSC.GT. 1 . )  COSSCM.O 
1  F( COSSC.LT.- 1 .  )  COSSC— 1.0 
SCANG=ACOS( COSSC ) 

SCATA=S1N(  .  i>*SCANO) 

l)E  NOM=  (  l .  +4 .  *  ( II  (  4 )  *  S  C  ATA )  * *  2  )  *  *  A  F.  X  P 

CoSX=SCALP( DCTO , SR ) 

XANG=ACOS( COSX ) 

S  INX=S  l.N(  XANG) 

CROSS=H(5)*(SINX**2)/ DENOM 

VOLTHT=CR0SS*EXPTG*EXPRG*EXP( -2 . *ATTEN) *WTHETA(  <) / ( ( SST*SSR) **2 ) 
C 

DOPPLER  SHIFT  CALCULATIONS 

DO  630  L=l ,3 
SD(L)-SR(L)-ST(L) 

VP( L)=VM( L) 

630  CONTINUE 

VP( 3)=H(6)+VP( 3) 

ANG=2 . *SCATA 
T1=H(7)*H(6)*H(6) 

V 1 =SCALP( VT , ST) 

V2=*SCALP(VR,SR) 

V3=SCALP(VP,SD) 

VS-V1-V2+V3 
DF*DFO 

DO  640  L-i.NFP 
PDQ=PDF(  DF , VS ,TI , ANG , BFREQ) 

V0LTH(  L)«V0LTI1(  L)+PDQ*VOLTHT 
DF=»DF+DELDF 
640  CONTINUE 

VOLZA-VOLZA+VOLTHT 

100  CONTINUE 

VOLRT=RRAD( J )*WRAD(  J ) 
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DO  650  L=1,NFP 

VOLR(L)=VOLR(L)+VOLTH( L)*VOLRT 
650  CONTINUE 

VOLZR=VOLZR+VOLZA*VOLRT 

200  CONTINUE 

VOLZT=PI*ZNORM*CONST*RPZ*RPZ*WZ(  I ) 

DO  660  L-l.NFP 

VOLZ( L)=VOLZ( L)+VOLR( L)*VOLZT 
660  CONTINUE 

VOLZZ=VOLZZ+VOLZR* VOLZT 

300  CONTINUE 

WRITE(6 , 730) 

730  FORMAT! 1H1.20X, 'VOLUME  INTEGRATED  SCATTERED  ENERGY  AT  DOPPLER- 
1SHIFTED  FREQUENCY '//4X, 'N' ,7X, 'DF,  KHZ' ,2X, ’DF*C/BFREQ' , 

2  2X, 'SCATTERED  POWER,  DB/HZ'/) 

DF=DF0 

DO  750  1=1, NFP 
FNORM=DF*C/ BFREQ 
DFKK.Z=DF/  1000. 

FVOLZ( I)=DFK4Z 
K=I+NFP 

IF( VOLZ( I) . LE.O. )  VOLZ( I )=1 .E-38 
F VOLZ ( K ) = 1 0 . *  ALOG 1 0 ( VOLZ ( I ) ) 

WRITE( 6, 740)  I , FV0LZ( I) ,FNORM,FVOLZ( K) 

DF=DF+DELDF 
750  CONTINUE 
VOLTDB=0. 

IF( VOLZZ . GT . 0 . )  VOLTDB* 1 0. *  ALOG 1 0 ( VOLZZ ) 

WRITE! 6, 760)  VOLTDB 

760  FORMAT! 1H  ,///40X, 'FREQUENCY  INTEGRATED  SCATTERED  POWER*  ’.F10.3, 
1  '  DB') 

740  FORMAT! 1H  ,2X, 13 , 2X, 3! 1PE1 1.4 , IX) ) 

PRINTER  PLOT  OF  DOPPLER  SHIFTED  SCATTERED  POWER 
CALL  PLOT! l,FVOLZ,NFP,2) 


RETURN 

END 

a***************************** ***** ** **** ****************************** 
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FUNCTION'  PDF(DF,  VMK , QQ ,  ANC ,  B FREQ) 

C 

GAUSSIAN  PROBABILITY  DENSITY  FUNCTION  FOR  DOPPLER  SHIFTED 
FREQUENCY 

AVERAGED  OVER  ALL  ANGLES 
DF=( F-BFREQ) 

VMK-EFFECTIVE  COMPONENT  OF  MEAN  RELATIVE  VELOVITY 
0/ “DIMENSIONAL  MEAN  SQUARE  TURBULENCE  INTENSITY 
ANC-MAGNITUDE  OF  PROJECTION  OF  SUM  OF  UNIT  VECTORS 

DATA  PI , C/3. 14 1592653, 2. 998E+08/ 


QOS-2 . *00* ANC* ANG 
PDF-0. 

IF!  OQS.EO.O.)  RETURN 
A»C*DF/ BFREQ-VMK 
ARG=A**2/QQS 
C0NST-SQRT(2.*PI*QQ) 

AVG-ABS(ANG) 

PDF»EXP( -ARG)*C/ ( CONST* AVG*BFREO) 

RETURN 

END 

FUNCTION  SCAT(AK) 

COMPUTES  0-FUNCTION  FOR  TOTAL  SCATTERING  CROSS-SECTION 


X- 2.*(A*K)**2 
A-CORRELATION  LENGTH 
K-WAVE  NUMBER 

C 

DATA  PI ,EX1 , EX2 , EX3/ 3. 1 4 1 592653 83333333 , .  1 6666667 , 1 .  1 666667/ 

SCAT-0. 

X«2.*(AK**2) 

XCUBE=X*<3 

IF( XCUBE. GT. 0, )  GO  TO  10 
RETURN 
10  CONTINUE 

XI- 1 .+2.*X 

SCAT=(2.*X*X+2.*X+1. )*( l.-Xl**EXl)/5.~ 

1  2. *( 1 .+X)*( X1**EX2- 1 . )+( X 1**EX3-1 . )/7 . 

SCAT=SCAT*6 . *P1/XCUBE 

RETURN 

END 

***★*★★*★** ***★**★*★**★*★*********★***********★*★*★*★****★*****★★****  *★ 
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SCALP-SCALAR  PRODUCT  OF  TWO  VECTORS  A  AND  B 

DIMENSION  A( 3) , B( 3) 

SCALP-0. 

DO  1  1=1,3 

SCALP- SCALP+A(  I)*B( I) 

1  CONTINUE 
RETURN 
END 

*******************************  **************************************** 


137 


NWC  TP  6386 


Appendix  C 

DERIVATION  OF  THE  RAY  TRACE  EQUATION 


This  appendix  outlines  the  derivation  of  the  equation  defining  the 
trajectory  of  the  refracted  rays  through  the  plume.  We  include  a 
derivation  that  starts  with  Maxwell's  equations  to  show  the 
approximations  made  in  the  geometric  optics  limit  of  EM  wave 
propagation  used  in  PRFIC. 

By  combining  Maxwell's  equation  for  the  electric  and  magnetic 
fields  t  and  A,  and  assuming  a  time  dependence  of  each  component  of  the 
form  e"*-u)t,  the  vector  wave  equation  (for  either  t  or  ft)  is  (Reference 
9) 


72£  +  k2n2£  -  -2Vt  •  -4  Vn2 

n* 


( C—  1 ) 


The  lowest  order  approximation  neglects  the  right  hand  side  for 
sufficiently  small  Vn2/n2.  The  approximate  equation  satisfied  by  the 
electric  field  is  then 


V2£  +  k2n2fi  -  0. 


(C-2) 


In  this  approximation,  the  components  of  the  electric  field  vector  all 
obey  the  same  equation,  and  it  is  possible  to  replace  the  vector  by  a 
scalar  amplitude,  $,  such  that 


V2$  +  k2n2$  *  0. 


(C-3) 


We  seek  solutions  to  this  equation  in  the  form 


<t>  -  AeikS,  (C-4) 


where  A,S  are  two  real  functions  of  position.  A  is  an  amplitude,  and  S 
is  a  phase,  called  the  Eikonel.  If  we  equate  real  and  imaginary  parts 
of  Eq .  (C-3)  for  $  given  by  Eq .  (C-4),  there  are  two  equations  for  A 
and  S 
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(VS)2  +  n2  -  0 


2V  —  +  V2S  -  0 
A 


(C-5) 


For  short  wavelengths  (large  k)  and  propagation  paths  such  that  V2 A/ A 
is  small,  the  first  equation  reduces  to  the  Eikonel  equation 


VS  -  n. 


(C-6) 


which  is  the  one  we  solve  to  determine  the  ray  trajectory.  The  second 
equation  gives  the  effect  of  ray  focusing  and  is  not  solved  in  this 
first  version  of  PRFIC.  It  reduces  to11 


1  Ji 

n2A  ds 


— — (n2A)  *  -V  •  t 


(C-7) 


and  its  solution  is 


n2A 


(n2A)Q 


V 


A 

t  ds 


(C-8) 


where  t  is  a  unit  vector  in  the  direction  of  the  ray  trajectory,  and  s 
is  the  length  along  the  trajectory.  A  is  essentially  the  area  of  a  ray 
bundle  propagating  along  the  path,  and  it  can  increase  _or  decrease 
according  to  the  integrated  divergence  of  the  unit  vector  t  along  the 
path.  Since  t  is  a  unit  vector,  it  has  no  divergence  except  that  due 
to  change  in  direction.  There  is  no  focusing  for  a  straight 
line-of-sight. 


^Molmud,  P. ,  "Raybend,  a  Ray  Tracing  Program  for  Microwave  Propagation 
Through  Rocket  Plumes,”  JANNAF  12c"  plume  Trajectory  Meeting, 
U.S.  Air  Force  Academy,  Colorado  Springs,  Nov.  18-20,  1980.  (paper, 

UNCLASSIFIED.)  also  (CPIA  Publication  332,  December  1980, 
pp.  273-312,  publication,  UNCLASSIFIED.) 
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Using  Che  geometry  of  the  wave  front  and  Its  normal,  the  Eikonel 
equation  can  be  manipulated  into  an  equation  for  the  ray  trajectory12 


( c-9  :• 


This  is  the  form  of  the  equation  solved  in  PRFIC.  The  solution 
proceeds  in  the  following  stepwise  manner  (Reference  7).  The  position 
vector  of  the  trajectory  is  t  -  jx(s),  y(s),  z(s)|  ,  and 


t 


r'(s) 


(C-10) 


The  solution  at  ?  (s  +  As)  is  obtained  from  that  at  t(s)  by  a 
series  expansion 

?(s  +  As)  -  r(s)  +  r'(s)As  +  y  r"(s)As2  ( C— 11) 

and 

r ' (s  +  as)  -  r'(s)  +  r"(s)  as  ( C— 1 2) 

By  geometry  we  obtain  (Reference  7) 

r"(s)  -  -  Vn  -  t(?  •  Vn)  (C-13) 

n 


Each  term  in  the  expansion  (Eq.  ( C— 11))  as  now  defined  at  s,  and  the 
values  at  s  +  As ,  can  be  obtained.  For  completeness,  Table  C-l 
provides  the  components  of  each  vector  in  the  rectangular  cartesian 
coordinate  system  used  in  the  calculation. 


12Kerr,  D.E.,  Propagation  of  Short  Radio  Waves,  McCraw  Hill,  New  York, 
1951,  pp.  41-44. 
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Appendix  D 

INDEX  OF  REFRACTION  FLUCTUATIONS 


The  purpose  of  this  appendix  is  to  outline  the  contribution  of  the 
fluctuating  plume  properties  to  the  index  of  refraction  structure 
constant 


Cn  "  1,6  JTT T  (D_1) 


is  the  ensemble  average  of  the  index  of  refraction  fluctuation. 
The  functional  dependence  of  n  is 

n  -  n  (up(ne),  ^(x^T)).  (D-2) 


A  Taylor  series  for  n*  can  be  written,  giving 
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In  principle,  the  mean  square  index  of  refraction  fluctuation  depends 
on  the  cross-correlation  of  the  electron  number  density,  temperature, 
and  chemical  species  fluctuations. 

This  degree  of  detail  is  not  generally  available  in  current  exhaust 
plume  flowfield  models.  The  g-equation  formulation  does  give  estimates 
of  these  quantities,  although  its  accuracy  in  afterburning,  ionized 
flows  has  not  been  established.  Moreover,  a  systematic  examination  of 
the  importance  of  each  term  in  Eq .  (D-4)  is  warranted  before  the  effort 
to  include  all  these  contributions  is  undertaken. 

In  the  verification  of  PRFIC,  we  have  used  only  the  first  term 
which  gives 


n  ’  2  -  n„ ' 2 
e 


>Y 


V  u>  / 

/ ven\2 " 

2 

4nJ 

1  +  - 

e 

\  w  / 

(D-5) 


or,  using  w2  -  nee2/e0mc,  oj/c  -  2tt/X  *  k,  and  the  Thompson  radius  of 
the  electron  rt  *  e2/Awe0mec;’ 
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n  “  TT 


"e"7  rl 
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V  \2 

en 


(D-6) 


We  have  used  an 
fluctuations  in  n, 


additional  assumption  for  this  study,  that  the 
,  depend  only  on  T  ,  through  the  equilibrium  (Saha) 


equation.  We  have  used  this  approach  for  the  purpose  of  verification, 
because  there  was  no  other  estimate  available  wi  thir^  the  SPF  code  at 
the  time.  We  do  have  a  g-equation  formulation  for  T  ,  based  on  the 


major  afterburning  species,  but  no  similar  results  for  n^ 
available  at  this  time. 


itself  were 


The  Saha  equation  for  the  electron  number  density  is 


ni  J2*mekTy/2 
n  \  h2  / 


exp  ~  W 


CD-7) 


(We  have  used  the  approximation  n  =  n  -  n_  in  the  denominator  because 
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r^/tx  «  1.)  With  this  relationship, 


and 


(D-8) 


(D-9) 


( D— 10) 


We  obtained  T*  from  our  SPF  modified  with  the  g-equation.  The 
ionization  potential,  I,  was  taken  to  be  that  of  potassium  (I/k  * 
50,354  K) .  The  actual  number  density  of  electrons,  ne ,  is  given  by  the 
kinetics  solution  in  SPF. 
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Appendix  E 

DOPPLER  SPECTRUM  OF  SCATTERED  POWER 


Radiation  scattered  by  the  turbulent  fluctuations  has  a  frequency 
spectrum  because  of  the  distribution  of  velocities  of  the  mean  and 
fluctuating  velocity.  The  frequency  shift  for  a  particular  velocity 
can  be  derived  by  a  sequential  application  of  the  Doppler  shift  to  the 
incident  and  scattered  wave.  At  the  scattering  volume,  the  apparent 
frequency  is  (see  Figure  E-l). 


v 


vo  ^T/S 


s 

c 


(E-l) 


and  at  the  receiver,  the  apparent  frequency  of  the  scattered  wave  is 


v 


V  ' 


+  V  ' 


(E-2) 


If  w  is  the  local  mean  axial  velocity  of  the  plume,  and  ?j'  the 
random  turbulent  fluctuation,  then 


-  iwk  +  $  +  q ' ) 


'T/S  =  vT  -  -T  vm  T  q 


(E-3) 


^S/R  *  <wk  +  Vra  +  q’>  -  trR 


(E-4) 


The  total  frequency  shift  between  the  transmitter  and  receiver  is 

then 


v"  =  v  +  —  ($T  •  s  -  $„  •  t)  +  (wk  +  $  )  •  ( t  -  s )  +  q ’  •  (t-s) 

( E-5) 
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or 


t)  +  (wk  +  0m) 


(s  -  O  +  q' 


(s  -  t) 
(E-6) 


Lot  the  component  of  the  fluctuating  velocity  that  contributes  a 
Doppler  shift  be  written  as 


(s  -  t) 


2u ' 


.  0 
sin  - 


( E— 7 ) 


where  u'  is  the  component  of  the  turbulent  fluctuating  velocity  in  the 
direction  (s  -  £).  At  each  point  in  the  plume,  there  will  be  a 
probability  distribution  of  frequencies  about  the  shift  due  to  the  mean 
velocities.  The  frequency  shift  due  only  to  the  mean  velocities  is 


*  ($x  *  s  -  •  t)  +  (wk  +  •  (s  -  t)  (E-8) 


For  the  first  version  of  PRF1C,  we  have  taken  the  probability 
density  function  of  the  velocity  fluctuations  to  be  Gaussian. 


P(u' ,  v' ,  w' )  = 


(2it  ^y'2 


exp  - 


u ' 2  +  v'2  +  w*  2 
2  q77 


(E-9) 


Furthermore,  we  assume  that  the  fluctuations  are  isotropic.  That  is, 
the  magnitude  of  the  total  fluctuating  velocity  is  independent  of 
direction.  With  the  assumption  of  isotropy,  we  can  choose  a  coordinate 
system  at  each  scattering  point  with  one  coordinate  direction  aligned 
with  the  direction  of  the  velocity  fluctuation  u'  in  Eq.  i'E-7). 
Fluctuations  in  the  other  two  directions  do  not  contribute  to  the 
Doppler  shift,  so  we  can  integrate  over  them  to  find  the  probability 
distribution  for  u' ,  alone. 
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P(u 


"-LL  P(u'  ,v' ,w' )dv'dw' 


(E-10) 


exp  - 


2q" 


(E-ll) 


With  the  relation  between  u*  and  Av  ( Eq .  (E-6)),  we  can  use  P(u') 
to  determine  P(Av) 


cAv 

v_ 


($T  •  s 


t)  -  (wk  +  $  )  •  (s  -  t) 

13 


(E-12) 


and 


using 


du’ 


cd( Av  ) 


2sin^  - 


(E-13) 


P(u')du'  “ ^  P(Av)dAv, 


(E-14) 


the  probability  density  for  the  Doppler  shifted  scattered  power  is 


P(  Av  ) 


•^(2"  q77) 


exp 


2v0  sin  - 


£^.vV 


8  q  '  2  sin*2  — 
2 


( E— 15) 
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with 


Vs  «  ($T  •  S  -  $R  •  t)  +  (wk  +  •  (s  -  t)  (E-16) 


the  effective  component  of  mean  velocities  contributing  to  a  Doppler 
shift.  All  the  individual  velocities  are  known  from  the  velocities  of 
the  transmitter,  receiver,  missile,  and  exhaust  plume  flowfield. 
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NOMENCLATURE 


A 

cn 

c 

E 

e 

Gt  *  Gr 

h 

I 

t. J.t 

K,  Kr,  Kj; 

Moo 

ne 

n 

ne 

P 

P(Av  ) 

P 

q 

Q 

r 

rt 

S 

5 

3 

T 

t 

u, v,w 
V 

x,  y,  z 
xi 


Coefficient  in  scalar  wave  equation 
Index  of  refraction  structure  constant 
Speed  of  light 
Electric  field  strength 
Electron  charge 

Transmitter  and  receiver  antenna  gains,  respectively 
Planck's  constant 
Ionization  potential 

Unit  vector  in  x,y,z  directions,  respectively 
Magnitude,  real,  and  imaginary  components 
of  complex  dielectric  constant 
Missile  flight  Mach  number 
Mass  of  electron 

Index  of  refraction,  or  total  particle  number  density 

Electron  number  density 

Power;  P0  is  incident  power 

Probability  density  distribution 

Pressure,  or  2k2  A2 

Root  mean  square  turbulent  velocity  fluctuation 

Total  scattering  function 

Radius,  or  range 

Thompson  electron  radius 

Eikonel 

Distance  along  ray  trajectory 

Unit  vector  from  scatterer  to  receiver 

Temperature 

Unit  vector  between  transmitter  and 
scattering  point 

Velocity  components  of  exhaust  plume 
Velocity  of  missile,  transmitter,  or  receiver 
Coordinate  directions 
Plume  species  mole  fractions 


a 


S 

"p 

10 


Attenuation  constant,  aspect  angle,  or 

direction  of  ray 

Phase  constant,  or  angle 

Angle 

Plasma  frequency,  rad/s 
Transmitter  frequency  rad/s 
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v 


X 


0 

4> 


X 

A 

e 


Frequency,  Hz  (v  =  io/2n) 

Collision  frequency 
Plasma  frequency,  Hz 
Permittivity  of  free  space 

Angle  between  electric  vector  and  scattering 

direction 

Scattering  angle 

Polar  coordinate 

Spectrum  function 

Scattering  cross  section 

Solid  angle 

Wavelength 

Turbulence  macroscale 
Turbulence  energy  dissipation 


Subscripts 

Transmitter,  receiver 
Chemical  species 
Electron 
Neutral 


T,  R 
i 
e 
n 


Superscripts 


Fluctuating  quantity 
Ensemble  average 
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